!!!! Business Cycles, Insurance - (10/22/2024) !!!!
!! Author: Michael Irwin
!! from "The Impact of Unemployment Insurance and Unsecured Credit on Business Cycles"
!! for Journal of Political Economy Macroeconomics
!! This code constains the replication to obtain the main results in the paper
!! Solution Method (individuals): OEGM for savings; grid search for debt
!! Solution Method (aggregate): Krusell and Smith 

!!!! State Variables !!!!
!! e : Idiosyncratic Productivity
!! a : Asset Level
!! j : Age
!! s : Credit Status
!! n : Individual Labor Market Status
!! b : Beta (discount factor)
!! z : Aggregate Exogenous State
    
!!!! Code Options !!!!
!! 1.) 'switch_scenario', determines which scenario is run
!! 2.) 'switch_results' collects results from simulation
    
!!!! Scenarios !!!!
!! 1.) Benchmark: extend duration of UI during recession
!! 2.) Decomposition: no fluctuations in TFP
!! 3.) Decomposition: no fluctuations in job separation rates
!! 4.) Decomposition: no fluctuations in job finding rates
!! 5.) Acyclical UI: no change in UI policies with the business cycle
!! 6.) UI Level Increase: counter-cyclical increase in replacement rate and max UI during recessions, budget neutral to UI extension
!! 7.) No Credit: model with no unsecured credit
!! 8.) No Credit, Acyclical UI
!! 9.) No Credit, UI Level Increase
!! 10.) Acyclical Credit: no aggregate fluctuations in terms of credit
!! 11.) Acyclical Credit, Acyclical UI
!! 12.) Acyclical Credit, UI Level Increase 
!! 13.) Recalibrate: no fluctuations in job finding rates

!! Code Introduction
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!! Model Parameters

program BCI_Main

use BCI_Par
use BCI_TB
use OEGM_TB
use Fortran_TB
use omp_lib

implicit none

!!!! Set the Scenario !!!!
integer, parameter :: switch_scenario=1

!!!! Variables for Forecasting Rules !!!!
integer :: count_g, count_b
real(rp) :: avg_K, sigma_K, avg_L, sigma_L, k_prime, l_prime, mean_K, mean_L
real(rp), dimension(nz) :: nu_k0, nu_k1, nu_k2, nu_l0, nu_l1, nu_l2
real(rp), dimension(nz) :: nu_p_k0, nu_p_k1, nu_p_k2, nu_p_l0, nu_p_l1, nu_p_l2
real(rp), dimension(nk) :: kgrid
real(rp), dimension(nl) :: lgrid
real(rp), dimension(3) :: b_g, b_b
integer, dimension(nm) :: z_index, k_index, l_index
real(rp), dimension(nt) :: L_sim, rand_sim
real(rp), dimension(nt+1) :: K_sim
real(rp), allocatable :: K_sim_g(:), K_sim_b(:), L_sim_g(:), L_sim_b(:), X_g(:,:), X_b(:,:)

!!!! Variables for Solving Model !!!!
integer :: iteration, count, seedsize, switch_results=0
integer :: a, e, j, n, s, b, i, p, t, a2, e2, j2, n2, i2, ia, ij, ibar_e, ibar_u, ip
integer :: z, z2, m, m2, k, k2, l, l2, iz, ik, il, im, im1, im2, im3, im4, iq
integer :: ngrids_ge, ngrids_gu, ngrids_be, ngrids_bu
real(rp) :: error, precision_outer
real(rp) :: rate, wage, rate_save
real(rp) :: Value, V_d, y_u, chi_e, chi_u, income_e, income_u
real(rp) :: omega, omega_k, omega_l
real(rp) :: maxV_e, maxG_e, maxC_e, maxD_e, maxQ_e, maxV_u, maxG_u, maxC_u, maxD_u, maxQ_u
real(rp) :: muval, gval, dval, hval, cval
real(rp) :: aggK, aggL, aggKp, avg_earn
real(rp) :: abar_e, abar_u

!!!! Variables for Results !!!!
integer, parameter :: n_income=nj*ni
integer, dimension(n_income) :: i_income, j_income
real(rp), dimension(n_income) :: income_grid, mu_income
real(rp) :: aggY, aggC, aggI, aggG, aggD, aggB, aggXi
real(rp) :: Dpop, Upop, Epop, Npop
real(rp) :: Dpop_E, Dpop_U, Dpop_P, Dpop_P_mean
real(rp) :: UI_exp, ret_exp, T_exp, tax_rev
real(rp) :: mean_Z, mean_lambda, mean_xi, mean_Y, mean_C, mean_I, mean_G, mean_D, mean_B, mean_U, mean_P, mean_r, mean_w, mean_M, mean_abar, mean_ebar
real(rp) :: sd_Z, sd_lambda, sd_xi, sd_Y, sd_C, sd_I, sd_G, sd_D, sd_B, sd_U, sd_P, sd_M, sd_abar
real(rp) :: corr_Y, corr_C, corr_I, corr_G, corr_D, corr_B, corr_U, corr_P, corr_M, corr_abar

!!!! Variables for Intermediary Profits !!!!
real(rp) :: qval, aggA, agg_discharge, profits

!!!! Variables for Borrowing Limits and Debt Prices !!!!
real(rp) :: avg_q, avg_limit, limit_good, pop_good, pop_borrow, blim

!!!! Variables for Individual Simulation !!!!
integer, parameter :: nsim=210
integer :: ih, ie, ix, it, j_return
integer, dimension(nj,nh,nsim) :: ie_sim, in_sim, ib_sim, ii_sim, ia_sim, is_sim
real(rp) :: omega_e, omega_a, omega_x
real(rp) :: c_dev, a_dev, mpc_mean, count_mpc
real(rp) :: earnings_minus, earnings_plus, debt_minus, debt_plus
real(rp) :: count_inc, count_dec, count_rr
real(rp) :: credit, credit_plus, credit_minus
real(rp) :: count_quarter, dc_mean, count_mean, return_i
real(rp), dimension(8) :: dc_q, dy_q, balance_q, bank_q
real(rp), dimension(ne) :: cum_pi_erg
real(rp), dimension(ne,ne) :: cum_pi_e
real(rp), allocatable :: rand_epsilon(:,:,:), rand_unemp(:,:,:), rand_beta(:,:,:), rand_credit(:,:,:), rand_upsilon(:,:,:), c_isim(:,:,:), h_isim(:,:,:), g_isim(:,:,:), d_isim(:,:,:), q_isim(:,:,:), inc_isim(:,:,:), abar_isim(:,:,:)

!!!! Arrays !!!!
integer, dimension(ni) :: n_index, e_index, b_index
integer, dimension(nm) :: Kp_index, Lp_index
integer, dimension(no) :: ilow_ge, ihigh_ge, ilow_gu, ihigh_gu, ilow_be, ihigh_be, ilow_bu, ihigh_bu
integer, dimension(nt) :: iz_sim
integer, dimension(ne,nn,nb) :: i_index
real(rp), dimension(na) :: agrid, EnGrid_ge, EnGrid_gu, EnGrid_be, EnGrid_bu
real(rp), dimension(na) :: EV_ge, EV_gu, EV_be, EV_bu, EVc_ge, EVc_gu, EVc_be, EVc_bu, DV_ge, DV_gu, DV_be, DV_bu
real(rp), dimension(nd) :: tempgrid, discount_e, discount_u
real(rp), dimension(ne) :: egrid, pi_erg
real(rp), dimension(nj) :: gamma_j
real(rp), dimension(nb) :: beta, phi_b
real(rp), dimension(nz) :: lambda, xi, zgrid, upsilon_d, upsilon_r, upsilon_bar, iota
real(rp), dimension(nm) :: r, w, E_r, E_rg, Kp_omega, Lp_omega
real(rp), dimension(nt) :: Z_sim, lambda_sim, xi_sim, Y_sim, C_sim, I_sim, G_sim, B_sim, D_sim, U_sim, P_sim, M_sim, r_sim, w_sim, abar_sim, ebar_sim
real(rp), dimension(nt-300) :: hp_Z, hp_lambda, hp_xi, hp_Y, hp_C, hp_I, hp_G, hp_D, hp_B, hp_U, hp_P, hp_M, hp_abar
real(rp), dimension(ne,ne) :: pi_eg, pi_eb
real(rp), dimension(nn,nn) :: pi_ng, pi_nb
real(rp), dimension(nz,nz) :: pi_z, cum_z
real(rp), dimension(nt-300,nt-300) :: HP, HP_inv
real(rp), dimension(nn,nn,nz) :: pi_n, cum_n
real(rp), dimension(ni,ni,2) :: pi_i, pi_ig, pi_ib
real(rp), allocatable :: mu(:,:,:,:), mu_next(:,:,:,:), mu_good(:,:,:,:), mu_bad(:,:,:,:), mur(:,:,:), mur_next(:,:,:)
real(rp), allocatable :: V(:,:,:,:,:), V_c(:,:,:,:,:), g(:,:,:,:,:), c(:,:,:,:,:), d(:,:,:,:,:), h(:,:,:,:,:), Vr(:,:,:,:), gr(:,:,:,:), cr(:,:,:,:), abar(:,:,:,:), q(:,:,:,:)
integer, allocatable :: seed(:)

!!!! Set the Number of Threads !!!!
!$ call omp_set_num_threads(np)

!!!! Variables for Aggregate Fluctuations !!!!
call agg_transitions(nz, rho_g, rho_b, pi_z, cum_z)
call aggregate_fluctuations(switch_scenario, zgrid, lambda, xi, iota)
call UI_policies(switch_scenario, upsilon_d, upsilon_r, upsilon_bar)

!!!! Forecasting Parameters for Capital !!!!
call forecasting_rules(switch_scenario, avg_K, sigma_K, nu_k0(:), nu_k1(:), nu_k2(:), avg_L, sigma_L, nu_L0(:), nu_L1(:), nu_L2(:))
    
!!!! Aggregate Indices !!!!
count = 0
do z = 1,nz
	do k = 1,nk
		do l = 1,nl
			count = count + 1
			z_index(count) = z
			k_index(count) = k
			l_index(count) = l
		end do
	end do
end do

!!!! Grid of Financial Wealth !!!!
call logspace(zero, -a_min, nd, tempgrid)
do a = 1,nd
	agrid(a) = -tempgrid(nd+1-a)
end do
call logspace(zero, a_max, na-nd+1, agrid(nd:na))

!!!! Persistent Earnings Process !!!!
call rouwenhorst(rho, sigma_eta_g, ne, pi_eg(:,:), egrid)
pi_eb(:,:) = pi_eg(:,:)
call ergodic(ne, pi_eg(:,:), pi_erg(:))
egrid = exp(egrid)
cum_pi_e(:,1) = pi_eg(:,1)
cum_pi_erg(1) = pi_erg(1)
do e = 2,ne
    cum_pi_e(:,e) = cum_pi_e(:,e-1) + pi_eg(:,e)
    cum_pi_erg(e) = cum_pi_erg(e-1) + pi_erg(e)
end do

!!!! Age-Component of Earnings !!!!
do j = 1,nj,4
	ij = ceiling(j/4.0_rp)
	gamma_j(j:j+3) = nu_1*dble(ij) + nu_2*dble(ij)**2.0_rp
end do
gamma_j = exp(gamma_j)

!!!! Transition of Employment States !!!!
call emp_transitions(nn, xi(1), lambda(1), upsilon_d(1), pi_ng(:,:))
call emp_transitions(nn, xi(2), lambda(2), upsilon_d(2), pi_nb(:,:))
pi_n(:,:,1) = pi_ng; pi_n(:,:,2) = pi_nb
cum_n(:,1,1) = pi_n(:,1,1); cum_n(:,1,2) = pi_n(:,1,2)
do n = 2,nn
    cum_n(:,n,1) = cum_n(:,n-1,1) + pi_n(:,n,1)
    cum_n(:,n,2) = cum_n(:,n-1,2) + pi_n(:,n,2)
end do

!!!! Heterogeneous Discount Factors !!!!
beta(1) = beta_h
beta(2) = beta_l
phi_b(1) = pi_beta
phi_b(2) = one-pi_beta

!!!! Assign Indices and Create Transition Matrix !!!!
call ind_transitions(ni, ne, nn, nb, pi_eg, pi_ng, e_index, n_index, b_index, i_index, pi_ig(:,:,:) )
call ind_transitions(ni, ne, nn, nb, pi_eb, pi_nb, e_index, n_index, b_index, i_index, pi_ib(:,:,:) )

!!!! Random Aggregate Shocks !!!!
seedsize = 0
call random_seed( size = seedsize )
allocate( seed(seedsize) )
seed = (/ 416641, 1112222 /)
call random_seed( put=seed )
call random_number(rand_sim)
iz_sim(1) = 1
do t = 2,nt
	iz = iz_sim(t-1)
	if ( rand_sim(t) <= cum_z(iz,1) ) then
		iz_sim(t) = 1
	else
		do z2 = 2,nz
			if ( rand_sim(t) > cum_z(iz,z2-1) .and. rand_sim(t) <= cum_z(iz,z2) ) then
				iz_sim(t) = z2
			end if
		end do	
	end if
end do

!!!! Initial Capital !!!!
K_sim(1) = avg_K

!!!! Set Income Distribution !!!!
call income_distribution(n_income, e_index, n_index, upsilon_r(1), gamma_j, egrid, i_income, j_income, income_grid)

!!Parameters of the Model
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!Begin Outer Loop

!!!! Outer Loop - Convergence of Forecasting Rules !!!!
error = 100.0_rp
iteration = 0

do while ( iteration <= max_iter )

    if ( iteration < max_iter ) then
        
        !!!! Set Grid for Aggregate Capital !!!!
        kgrid(1) = avg_K - 3.0_rp*sigma_K
        kgrid(nk) = avg_K + 3.0_rp*sigma_K
        call linspace(kgrid(1), kgrid(nk), nk, kgrid(:))

        !!!! Set Grid for Aggregate Labor !!!!
        lgrid(1) = avg_L - 3.0_rp*sigma_L
        lgrid(nl) = avg_L + 3.0_rp*sigma_L
        call linspace(lgrid(1), lgrid(nl), nl, lgrid(:))
        
    end if

    !!!! Functions Dependent on Aggregate State Space !!!!
    do m = 1,nm
	    z = z_index(m); k = k_index(m); l = l_index(m)
        r(m) = zgrid(z)*alpha*( kgrid(k)**(alpha-one) )*( lgrid(l)**(one-alpha) ) - delta
        w(m) = zgrid(z)*(one-alpha)*( kgrid(k)**alpha )*( lgrid(l)**(-alpha) )
    
	    k_prime = exp( nu_k0(z) + nu_k1(z)*log(kgrid(k)) + nu_k2(z)*log(lgrid(l)) )
	    call interpolation( kgrid, k_prime, 1, nk, nk, ik, omega_k )
        Kp_index(m) = ik
        Kp_omega(m) = omega_k
    
	    l_prime = exp( nu_l0(z) + nu_l1(z)*log(lgrid(l)) + nu_l2(z)*log(kgrid(k)) )
	    call interpolation( lgrid, l_prime, 1, nl, nl, il, omega_l )
        Lp_index(m) = il
        Lp_omega(m) = omega_l
    
        E_r(m) = pi_z(z,1)*( zgrid(1)*alpha*( k_prime**(alpha-one) )*( l_prime**(one-alpha) ) - delta ) + pi_z(z,2)*( zgrid(2)*alpha*( k_prime**(alpha-one) )*( l_prime**(one-alpha) ) - delta )
        if ( z == 1 ) then
            E_rg(m) = E_r(m)
        else
            E_rg(m) = E_r(m-49)
        end if
    end do

!!Begin Outer Loop
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!Solving for Decision Rules 

    !!!! Allocate Large Arrays !!!!
    if ( iteration < 1 ) then
	    allocate(V(ni,na,nj,ns,nm))
	    allocate(g(ni,na,nj,ns,nm))
	    allocate(c(ni,na,nj,ns,nm))
	    allocate(d(ni,na,nj,ns,nm))
	    allocate(h(ni,na,nj,ns,nm))
	    allocate(Vr(nb,na,nR,nm))
	    allocate(gr(nb,na,nR,nm))
	    allocate(cr(nb,na,nR,nm))
	    allocate(abar(ni,na,nj,nm))
	    allocate(q(ni,na,nj,nm))
        allocate(V_c(ni,na,nj,ns,nm))
    end if

    !!!! Initial Functional Values !!!!
    V = zero; g=zero; c=zero; d=zero; h=zero; Vr=zero; gr=zero; cr=zero
    abar=zero; q=zero
    V_c = zero

    !!!! Final Age Value !!!!
    do a = nd,na
        Vr(:,a,nR,:) = (one/(one-sigma))*(( agrid(a) + y_r )**(one-sigma))
        cr(:,a,nR,:) = agrid(a) + y_r
        gr(:,a,nR,:) = zero
    end do

    !!!! Decision Rules in Retirement - EGM !!!!
    do j = (nR-1),1,-1
	    do b = 1,nb
		    do m = 1,nm
			    z = z_index(m); k = k_index(m); l = l_index(m)
                rate = r(m); wage = w(m); rate_save = E_r(m)
                ik = Kp_index(m); omega_k = Kp_omega(m); il = Lp_index(m); omega_l = Lp_omega(m)

			    !!!! Derivatives of Expected Future Value !!!!
			    EV_gu = zero
			    do a2 = na,nd,-1
				    Value = zero
				    do m2 = 1,nm
					    z2 = z_index(m2); k2 = k_index(m2); l2 = l_index(m2)

					    if ( k2 == ik .and. l2 == il ) then
						    Value = Value + omega_k*omega_l*pi_z(z,z2)*Vr(b,a2,j+1,m2)
					    elseif ( k2 == ik+1 .and. l2 == il ) then
						    Value = Value + (one-omega_k)*omega_l*pi_z(z,z2)*Vr(b,a2,j+1,m2)
					    elseif ( k2 == ik .and. l2 == il+1 ) then
						    Value = Value + omega_k*(one-omega_l)*pi_z(z,z2)*Vr(b,a2,j+1,m2)
					    elseif ( k2 == ik+1 .and. l2 == il+1 ) then
						    Value = Value + (one-omega_k)*(one-omega_l)*pi_z(z,z2)*Vr(b,a2,j+1,m2)
					    end if
				    end do
				    EV_gu(a2) = Value

				    if ( a2 < na ) then
					    DV_gu(a2) = ( EV_gu(a2+1)-EV_gu(a2) )/( agrid(a2+1)-agrid(a2) )
                    end if
			    end do
			    DV_gu(na) = DV_gu(na-1)

			    !!!! Solve for Decision Rules !!!!
                call endogenous_grids(beta(b), rate_save, y_r - y_bar, agrid, DV_gu, ngrids_gu, ilow_gu, ihigh_gu, EnGrid_gu)
                maxV_u = disutility; maxG_u = zero; maxC_u = zero; maxQ_u = zero
                do a = nd,na
				    call oegm_interpolation(ngrids_gu, beta(b), zero, agrid(a), y_r - y_bar, rate_save, ilow_gu, ihigh_gu, EnGrid_gu, agrid, EV_gu, maxV_u, maxG_u, maxC_u, maxQ_u)
				    Vr(b,a,j,m) = maxV_u
				    gr(b,a,j,m) = maxG_u
				    cr(b,a,j,m) = maxC_u
                end do

		    end do
	    end do
    end do

    !!!! Decision Rules Working Ages - OEGM !!!!
    do j = nj,1,-1

	    !$omp parallel do private(i, z, z2, k, k2, l, l2, m, m2, a, a2, n, e, b, s, discount_e, discount_u, ibar_e, ibar_u, abar_e, abar_u, pi_i, &
        !$              EV_ge, EV_gu, EV_be, EV_bu, EVc_ge, EVc_gu, EVc_be, EVc_bu, DV_ge, DV_gu, DV_be, DV_bu, EnGrid_ge, EnGrid_gu, EnGrid_be, EnGrid_bu, &
        !$              ilow_ge, ilow_gu, ilow_be, ilow_bu, ihigh_ge, ihigh_gu, ihigh_be, ihigh_bu, ngrids_ge, ngrids_gu, ngrids_be, ngrids_bu, &
	    !$              maxV_e, maxG_e, maxC_e, maxD_e, maxQ_e, maxV_u, maxG_u, maxC_u, maxD_u, maxQ_u, p, ia, omega, &
	    !$              Value, V_d, y_u, ik, il, omega_k, omega_l, k_prime, l_prime, rate, wage, rate_save, income_e, income_u, chi_e, chi_u )
	
	    do i = 1,ni
		    n = n_index(i); e = e_index(i); b = b_index(i) 
        
		    do m = 1,nm
                z = z_index(m); k = k_index(m); l = l_index(m)
                rate = r(m); wage = w(m); rate_save = E_r(m)
                ik = Kp_index(m); omega_k = Kp_omega(m); il = Lp_index(m); omega_l = Lp_omega(m)
            
                !!!! Variables Dependent on Labor Choice !!!!
                if ( n == 1 ) then
                    income_e = (one-tau)*wage*egrid(e)*gamma_j(j)
                    income_u = (one-tau)*wage*egrid(e)*gamma_j(j)
                    chi_e = chi_w
                    chi_u = chi_w
                elseif ( n == 2 ) then
                    if ( upsilon_r(z)*wage*egrid(e)*gamma_j(j) >= upsilon_bar(z) ) then
                        income_e = upsilon_bar(z)
                    else
            		    income_e = upsilon_r(z)*wage*egrid(e)*gamma_j(j)
                    end if
                    income_u = zero
                    chi_e = chi_s
                    chi_u = zero
                else
                    income_e = zero
                    income_u = zero
                    chi_e = chi_s
                    chi_u = zero
                end if

			    if ( z == 1 ) then
                    pi_i = pi_ig
                else
                    pi_i = pi_ib
                end if
            
                !!!! Calculate Expected Future Variables !!!!
                if ( j == nj ) then
                
                    discount_e(:) = zero; discount_u(:) = zero
                    ibar_e = nd; ibar_u = nd
                    abar_e = zero; abar_u = zero
                    EV_ge = zero; EV_gu = zero; EV_be = zero; EV_bu = zero
                    EVc_ge = zero; EVc_gu = zero; EVc_be = zero; EVc_bu = zero
                
    		        do a2 = na,nd,-1
                        Value = zero
                        do m2 = 1,nm
                	        z2 = z_index(m2); k2 = k_index(m2); l2 = l_index(m2)
                            if ( k2 == ik .and. l2 == il ) then
                        	    Value = Value + omega_k*omega_l*pi_z(z,z2)*Vr(b,a2,1,m2)
                            elseif ( k2 == ik+1 .and. l2 == il ) then
                        	    Value = Value + (one-omega_k)*omega_l*pi_z(z,z2)*Vr(b,a2,1,m2)
                            elseif ( k2 == ik .and. l2 == il+1 ) then
                            	Value = Value + omega_k*(one-omega_l)*pi_z(z,z2)*Vr(b,a2,1,m2)
                            elseif ( k2 == ik+1 .and. l2 == il+1 ) then
                        	    Value = Value + (one-omega_k)*(one-omega_l)*pi_z(z,z2)*Vr(b,a2,1,m2)
                            end if
                        end do
                	    EV_ge(a2) = Value; EV_gu(a2) = Value; EV_be(a2) = Value; EV_bu(a2) = Value
                        EVc_ge(a2) = Value; EVc_gu(a2) = Value; EVc_be(a2) = Value; EVc_bu(a2) = Value

                        if ( a2 < na ) then
                    	    DV_ge(a2) = ( EV_ge(a2+1)-EV_ge(a2) )/( agrid(a2+1)-agrid(a2) )
                	        DV_gu(a2) = ( EV_gu(a2+1)-EV_gu(a2) )/( agrid(a2+1)-agrid(a2) )
                	        DV_be(a2) = ( EV_be(a2+1)-EV_be(a2) )/( agrid(a2+1)-agrid(a2) )
                    	    DV_bu(a2) = ( EV_bu(a2+1)-EV_bu(a2) )/( agrid(a2+1)-agrid(a2) )
                        end if
                    
        	        end do
                    DV_ge(na) = DV_ge(na-1); DV_gu(na) = DV_gu(na-1); DV_be(na) = DV_be(na-1); DV_bu(na) = DV_bu(na-1)
                
                else
                
                    !!!! Discount Prices !!!!
                    if ( switch_scenario >= 10 .and. switch_scenario <= 12 ) then           !Acyclial Credit Scenarios
                        discount_e(:) = (one-iota(1))/(one+E_rg(m))
                        discount_u(:) = discount_e(:)
                        do a2 = 1,(nd-1)
				            do m2 = 1,nm
					            z2 = z_index(m2); k2 = k_index(m2); l2 = l_index(m2)
					            if ( k2 == ik .and. l2 == il ) then
						            discount_e(a2) = discount_e(a2) - omega_k*omega_l*pi_z(1,z2)*dot_product(pi_ig(i,:,1),d(:,a2,j+1,1,m2))/(one+E_rg(m))
    						        discount_u(a2) = discount_u(a2) - omega_k*omega_l*pi_z(1,z2)*dot_product(pi_ig(i,:,2),d(:,a2,j+1,1,m2))/(one+E_rg(m))
	    				        elseif ( k2 == ik+1 .and. l2 == il ) then
		    				        discount_e(a2) = discount_e(a2) - (one-omega_k)*omega_l*pi_z(1,z2)*dot_product(pi_ig(i,:,1),d(:,a2,j+1,1,m2))/(one+E_rg(m))
			    			        discount_u(a2) = discount_u(a2) - (one-omega_k)*omega_l*pi_z(1,z2)*dot_product(pi_ig(i,:,2),d(:,a2,j+1,1,m2))/(one+E_rg(m))
				    	        elseif ( k2 == ik .and. l2 == il+1 ) then
					    	        discount_e(a2) = discount_e(a2) - omega_k*(one-omega_l)*pi_z(1,z2)*dot_product(pi_ig(i,:,1),d(:,a2,j+1,1,m2))/(one+E_rg(m))
						            discount_u(a2) = discount_u(a2) - omega_k*(one-omega_l)*pi_z(1,z2)*dot_product(pi_ig(i,:,2),d(:,a2,j+1,1,m2))/(one+E_rg(m))
					            elseif ( k2 == ik+1 .and. l2 == il+1 ) then
						            discount_e(a2) = discount_e(a2) - (one-omega_k)*(one-omega_l)*pi_z(1,z2)*dot_product(pi_ig(i,:,1),d(:,a2,j+1,1,m2))/(one+E_rg(m))
    						        discount_u(a2) = discount_u(a2) - (one-omega_k)*(one-omega_l)*pi_z(1,z2)*dot_product(pi_ig(i,:,2),d(:,a2,j+1,1,m2))/(one+E_rg(m))
	    				        end if
                           end do
                        end do
                    else
                        discount_e(:) = (one-iota(z))/(one+rate_save)
                        discount_u(:) = discount_e(:)
	    		        do a2 = 1,(nd-1)
		    		        do m2 = 1,nm
			    		        z2 = z_index(m2); k2 = k_index(m2); l2 = l_index(m2)
				    	        if ( k2 == ik .and. l2 == il ) then
					    	        discount_e(a2) = discount_e(a2) - omega_k*omega_l*pi_z(z,z2)*dot_product(pi_i(i,:,1),d(:,a2,j+1,1,m2))/(one+rate_save)
						            discount_u(a2) = discount_u(a2) - omega_k*omega_l*pi_z(z,z2)*dot_product(pi_i(i,:,2),d(:,a2,j+1,1,m2))/(one+rate_save)
    					        elseif ( k2 == ik+1 .and. l2 == il ) then
	    					        discount_e(a2) = discount_e(a2) - (one-omega_k)*omega_l*pi_z(z,z2)*dot_product(pi_i(i,:,1),d(:,a2,j+1,1,m2))/(one+rate_save)
		    				        discount_u(a2) = discount_u(a2) - (one-omega_k)*omega_l*pi_z(z,z2)*dot_product(pi_i(i,:,2),d(:,a2,j+1,1,m2))/(one+rate_save)
			    		        elseif ( k2 == ik .and. l2 == il+1 ) then
				    		        discount_e(a2) = discount_e(a2) - omega_k*(one-omega_l)*pi_z(z,z2)*dot_product(pi_i(i,:,1),d(:,a2,j+1,1,m2))/(one+rate_save)
					    	        discount_u(a2) = discount_u(a2) - omega_k*(one-omega_l)*pi_z(z,z2)*dot_product(pi_i(i,:,2),d(:,a2,j+1,1,m2))/(one+rate_save)
					            elseif ( k2 == ik+1 .and. l2 == il+1 ) then
						            discount_e(a2) = discount_e(a2) - (one-omega_k)*(one-omega_l)*pi_z(z,z2)*dot_product(pi_i(i,:,1),d(:,a2,j+1,1,m2))/(one+rate_save)
						            discount_u(a2) = discount_u(a2) - (one-omega_k)*(one-omega_l)*pi_z(z,z2)*dot_product(pi_i(i,:,2),d(:,a2,j+1,1,m2))/(one+rate_save)
			    		        end if
		    		        end do
                        end do
                    end if

    			    !!!! Calculate the Endogenous Borrowing Limit !!!!
                    if ( switch_scenario == 7 .or. switch_scenario == 8 .or. switch_scenario == 9 .or. switch_scenario == 15 ) then         !No Credit Scenarios
                        ibar_e = nd
                        ibar_u = nd
                    else
			            call borrowing_limit(nd, agrid(1:nd), discount_e(:), ibar_e, abar_e)
    			        call borrowing_limit(nd, agrid(1:nd), discount_u(:), ibar_u, abar_u)
                    end if
            
                    !!!! Calculate the Expected Future Value and Numerical Derivatives !!!!
                    EV_ge = zero; EV_gu = zero; EV_be = zero; EV_bu = zero
                    do a2 = na,1,-1	
                	    do m2 = 1,nm
                		    z2 = z_index(m2); k2 = k_index(m2); l2 = l_index(m2)
            	    	    if ( k2 == ik .and. l2 == il ) then
                	    	    EV_ge(a2) = EV_ge(a2) + omega_k*omega_l*pi_z(z,z2)*dot_product(pi_i(i,:,1),V(:,a2,j+1,1,m2))
                    		    EV_gu(a2) = EV_gu(a2) + omega_k*omega_l*pi_z(z,z2)*dot_product(pi_i(i,:,2),V(:,a2,j+1,1,m2))
                    		    EV_be(a2) = EV_be(a2) + omega_k*omega_l*pi_z(z,z2)*dot_product(pi_i(i,:,1),V(:,a2,j+1,2,m2))
                    		    EV_bu(a2) = EV_bu(a2) + omega_k*omega_l*pi_z(z,z2)*dot_product(pi_i(i,:,2),V(:,a2,j+1,2,m2))
                    	    elseif ( k2 == ik+1 .and. l2 == il ) then
                    		    EV_ge(a2) = EV_ge(a2) + (one-omega_k)*omega_l*pi_z(z,z2)*dot_product(pi_i(i,:,1),V(:,a2,j+1,1,m2))
                	    	    EV_gu(a2) = EV_gu(a2) + (one-omega_k)*omega_l*pi_z(z,z2)*dot_product(pi_i(i,:,2),V(:,a2,j+1,1,m2))
                		        EV_be(a2) = EV_be(a2) + (one-omega_k)*omega_l*pi_z(z,z2)*dot_product(pi_i(i,:,1),V(:,a2,j+1,2,m2))
                		        EV_bu(a2) = EV_bu(a2) + (one-omega_k)*omega_l*pi_z(z,z2)*dot_product(pi_i(i,:,2),V(:,a2,j+1,2,m2))
                    	    elseif ( k2 == ik .and. l2 == il+1 ) then
                    		    EV_ge(a2) = EV_ge(a2) + omega_k*(one-omega_l)*pi_z(z,z2)*dot_product(pi_i(i,:,1),V(:,a2,j+1,1,m2))
                    		    EV_gu(a2) = EV_gu(a2) + omega_k*(one-omega_l)*pi_z(z,z2)*dot_product(pi_i(i,:,2),V(:,a2,j+1,1,m2))
                    		    EV_be(a2) = EV_be(a2) + omega_k*(one-omega_l)*pi_z(z,z2)*dot_product(pi_i(i,:,1),V(:,a2,j+1,2,m2))
                    		    EV_bu(a2) = EV_bu(a2) + omega_k*(one-omega_l)*pi_z(z,z2)*dot_product(pi_i(i,:,2),V(:,a2,j+1,2,m2))
                	        elseif ( k2 == ik+1 .and. l2 == il+1 ) then
                		        EV_ge(a2) = EV_ge(a2) + (one-omega_k)*(one-omega_l)*pi_z(z,z2)*dot_product(pi_i(i,:,1),V(:,a2,j+1,1,m2))
                    		    EV_gu(a2) = EV_gu(a2) + (one-omega_k)*(one-omega_l)*pi_z(z,z2)*dot_product(pi_i(i,:,2),V(:,a2,j+1,1,m2))
                    		    EV_be(a2) = EV_be(a2) + (one-omega_k)*(one-omega_l)*pi_z(z,z2)*dot_product(pi_i(i,:,1),V(:,a2,j+1,2,m2))
                    		    EV_bu(a2) = EV_bu(a2) + (one-omega_k)*(one-omega_l)*pi_z(z,z2)*dot_product(pi_i(i,:,2),V(:,a2,j+1,2,m2))
                    	    end if
                        end do	

            		    if ( a2 < na ) then
            			    DV_ge(a2) = ( EV_ge(a2+1)-EV_ge(a2) )/( agrid(a2+1)-agrid(a2) )
            			    DV_gu(a2) = ( EV_gu(a2+1)-EV_gu(a2) )/( agrid(a2+1)-agrid(a2) )
        	    		    DV_be(a2) = ( EV_be(a2+1)-EV_be(a2) )/( agrid(a2+1)-agrid(a2) )
        		    	    DV_bu(a2) = ( EV_bu(a2+1)-EV_bu(a2) )/( agrid(a2+1)-agrid(a2) )
        		        end if
            	    end do
            	    DV_ge(na) = DV_ge(na-1); DV_gu(na) = DV_gu(na-1); DV_be(na) = DV_be(na-1); DV_bu(na) = DV_bu(na-1)
            
                    !!!! Expected Value of Consumption - Welfare Analysis !!!!
                    EVc_ge = zero; EVc_gu = zero; EVc_be = zero; EVc_bu = zero
                    do a2 = na,1,-1	
        	            do m2 = 1,nm
        		            z2 = z_index(m2); k2 = k_index(m2); l2 = l_index(m2)
        		            if ( k2 == ik .and. l2 == il ) then
                                EVc_ge(a2) = EVc_ge(a2) + omega_k*omega_l*pi_z(z,z2)*dot_product(pi_i(i,:,1),V_c(:,a2,j+1,1,m2))
                                EVc_gu(a2) = EVc_gu(a2) + omega_k*omega_l*pi_z(z,z2)*dot_product(pi_i(i,:,2),V_c(:,a2,j+1,1,m2))
                                EVc_be(a2) = EVc_be(a2) + omega_k*omega_l*pi_z(z,z2)*dot_product(pi_i(i,:,1),V_c(:,a2,j+1,2,m2))
                    	        EVc_bu(a2) = EVc_bu(a2) + omega_k*omega_l*pi_z(z,z2)*dot_product(pi_i(i,:,2),V_c(:,a2,j+1,2,m2))
                	        elseif ( k2 == ik+1 .and. l2 == il ) then
                		        EVc_ge(a2) = EVc_ge(a2) + (one-omega_k)*omega_l*pi_z(z,z2)*dot_product(pi_i(i,:,1),V_c(:,a2,j+1,1,m2))
            	    	        EVc_gu(a2) = EVc_gu(a2) + (one-omega_k)*omega_l*pi_z(z,z2)*dot_product(pi_i(i,:,2),V_c(:,a2,j+1,1,m2))
                                EVc_be(a2) = EVc_be(a2) + (one-omega_k)*omega_l*pi_z(z,z2)*dot_product(pi_i(i,:,1),V_c(:,a2,j+1,2,m2))
                                EVc_bu(a2) = EVc_bu(a2) + (one-omega_k)*omega_l*pi_z(z,z2)*dot_product(pi_i(i,:,2),V_c(:,a2,j+1,2,m2))
                            elseif ( k2 == ik .and. l2 == il+1 ) then
                    	        EVc_ge(a2) = EVc_ge(a2) + omega_k*(one-omega_l)*pi_z(z,z2)*dot_product(pi_i(i,:,1),V_c(:,a2,j+1,1,m2))
                		        EVc_gu(a2) = EVc_gu(a2) + omega_k*(one-omega_l)*pi_z(z,z2)*dot_product(pi_i(i,:,2),V_c(:,a2,j+1,1,m2))
                		        EVc_be(a2) = EVc_be(a2) + omega_k*(one-omega_l)*pi_z(z,z2)*dot_product(pi_i(i,:,1),V_c(:,a2,j+1,2,m2))
            	    	        EVc_bu(a2) = EVc_bu(a2) + omega_k*(one-omega_l)*pi_z(z,z2)*dot_product(pi_i(i,:,2),V_c(:,a2,j+1,2,m2))
            	            elseif ( k2 == ik+1 .and. l2 == il+1 ) then
                                EVc_ge(a2) = EVc_ge(a2) + (one-omega_k)*(one-omega_l)*pi_z(z,z2)*dot_product(pi_i(i,:,1),V_c(:,a2,j+1,1,m2))
                                EVc_gu(a2) = EVc_gu(a2) + (one-omega_k)*(one-omega_l)*pi_z(z,z2)*dot_product(pi_i(i,:,2),V_c(:,a2,j+1,1,m2))
                                EVc_be(a2) = EVc_be(a2) + (one-omega_k)*(one-omega_l)*pi_z(z,z2)*dot_product(pi_i(i,:,1),V_c(:,a2,j+1,2,m2))
                    	        EVc_bu(a2) = EVc_bu(a2) + (one-omega_k)*(one-omega_l)*pi_z(z,z2)*dot_product(pi_i(i,:,2),V_c(:,a2,j+1,2,m2))
                	        end if
                        end do	
                    end do
            
                end if

                !!!! Form Endogenous Grids !!!!
                call endogenous_grids(beta(b), rate_save, income_e, agrid, DV_ge, ngrids_ge, ilow_ge, ihigh_ge, EnGrid_ge)
                call endogenous_grids(beta(b), rate_save, income_u, agrid, DV_gu, ngrids_gu, ilow_gu, ihigh_gu, EnGrid_gu)
                call endogenous_grids(beta(b), rate_save, income_e, agrid, theta*DV_ge(:) + (one-theta)*DV_be(:), ngrids_be, ilow_be, ihigh_be, EnGrid_be)
                call endogenous_grids(beta(b), rate_save, income_u, agrid, theta*DV_gu(:) + (one-theta)*DV_bu(:), ngrids_bu, ilow_bu, ihigh_bu, EnGrid_bu)

                !!!! Solve for Decision Rules !!!!
                do a = 1,na

                    !!!! Good Credit - Work/Search !!!!
                    s = 1
                    call grid_search(ibar_e, beta(b), chi_e, agrid(a), income_e, discount_e, agrid, EV_ge, maxV_e, maxG_e, maxC_e, maxQ_e)
                    call oegm_interpolation(ngrids_ge, beta(b), chi_e, agrid(a), income_e, rate_save, ilow_ge, ihigh_ge, EnGrid_ge, agrid, EV_ge, maxV_e, maxG_e, maxC_e, maxQ_e)
                    V_d = (one/(one-sigma))*(( income_e + y_bar )**(one-sigma)) - chi_e - chi_b + beta(b)*EV_be(nd)

                    !!!! Bankruptcy Decision - Work/Search !!!!
                    if ( V_d > maxV_e ) then   
                    	maxV_e = V_d
                    	maxG_e = zero
                    	maxC_e = income_e + y_bar
                    	maxD_e = one
                        maxQ_e = zero
                    else
                    	maxD_e = zero
                    end if

                    !!!! Good Credit - Leisure !!!!
                    call grid_search(ibar_u, beta(b), chi_u, agrid(a), income_u, discount_u, agrid, EV_gu, maxV_u, maxG_u, maxC_u, maxQ_u)
                    call oegm_interpolation(ngrids_gu, beta(b), chi_u, agrid(a), income_u, rate_save, ilow_gu, ihigh_gu, EnGrid_gu, agrid, EV_gu, maxV_u, maxG_u, maxC_u, maxQ_u)
                    V_d = (one/(one-sigma))*(( income_u + y_bar )**(one-sigma)) - chi_u - chi_b + beta(b)*EV_bu(nd)

                    !!!! Bankruptcy Decision - Leisure !!!!
                    if ( V_d > maxV_u ) then
                    	maxV_u = V_d
                    	maxG_u = zero
                    	maxC_u = income_u + y_bar
                	    maxD_u = one
                        maxQ_u = zero
                    else
                        maxD_u = zero
                    end if

                    !!!! Labor Market Decision !!!!
                	if ( maxV_e >= maxV_u ) then
            			V(i,a,j,s,m) = maxV_e
    		    		g(i,a,j,s,m) = maxG_e
                        c(i,a,j,s,m) = maxC_e
                	    d(i,a,j,s,m) = maxD_e
                		h(i,a,j,s,m) = one
                        q(i,a,j,m) = maxQ_e
                        abar(i,a,j,m) = abar_e
                        call interpolation(agrid, maxG_e, 1, na, na, ia, omega)
                        V_c(i,a,j,s,m) = (one/(one-sigma))*(( maxC_e )**(one-sigma)) + beta(b)*( omega*EVc_ge(ia) + (one-omega)*EVc_ge(ia+1) )    
				    else
                        V(i,a,j,s,m) = maxV_u
                    	g(i,a,j,s,m) = maxG_u
                		c(i,a,j,s,m) = maxC_u
        	    		d(i,a,j,s,m) = maxD_u
    			    	h(i,a,j,s,m) = zero
                        q(i,a,j,m) = maxQ_u
                        abar(i,a,j,m) = abar_u
                        call interpolation(agrid, maxG_u, 1, na, na, ia, omega)
                        V_c(i,a,j,s,m) = (one/(one-sigma))*(( maxC_u )**(one-sigma)) + beta(b)*( omega*EVc_gu(ia) + (one-omega)*EVc_gu(ia+1) )
                        
        		    end if

                    !!!! Bad Credit Status !!!!
                    s = 2
                    if ( a >= nd )  then

                        !!!! Decision Rules from Working !!!!
                    	maxV_e = disutility; maxG_e = zero; maxC_e = zero; maxQ_e = zero
                		call oegm_interpolation(ngrids_be, beta(b), chi_e, agrid(a), income_e, rate_save, ilow_be, ihigh_be, EnGrid_be, agrid, theta*EV_ge(:)+(one-theta)*EV_be(:), maxV_e, maxG_e, maxC_e, maxQ_e)

                        !!!! Start Decision Rules from Leisure !!!!
                    	maxV_u = disutility; maxG_u = zero; maxC_u = zero; maxQ_u = zero
                		call oegm_interpolation(ngrids_bu, beta(b), chi_u, agrid(a), income_u, rate_save, ilow_bu, ihigh_bu, EnGrid_bu, agrid, theta*EV_gu(:)+(one-theta)*EV_bu(:), maxV_u, maxG_u, maxC_u, maxQ_u)
    
                        !!!! Labor Supply Decision !!!!
                    	if ( maxV_e >= maxV_u ) then
                			V(i,a,j,s,m) = maxV_e
        		    		g(i,a,j,s,m) = maxG_e
                        	c(i,a,j,s,m) = maxC_e
                    	    d(i,a,j,s,m) = zero
                        	h(i,a,j,s,m) = one
                            call interpolation(agrid, maxG_e, 1, na, na, ia, omega)
                            V_c(i,a,j,s,m) = (one/(one-sigma))*(( maxC_e )**(one-sigma)) + beta(b)*theta*( omega*EVc_ge(ia) + (one-omega)*EVc_ge(ia+1) ) + beta(b)*(one-theta)*( omega*EVc_be(ia) + (one-omega)*EVc_be(ia+1) )            
                        else
                        	V(i,a,j,s,m) = maxV_u
                        	g(i,a,j,s,m) = maxG_u
                        	c(i,a,j,s,m) = maxC_u
                        	d(i,a,j,s,m) = zero
                        	h(i,a,j,s,m) = zero
                            call interpolation(agrid, maxG_u, 1, na, na, ia, omega)
                            V_c(i,a,j,s,m) = (one/(one-sigma))*(( maxC_u )**(one-sigma)) + beta(b)*theta*( omega*EVc_gu(ia) + (one-omega)*EVc_gu(ia+1) ) + beta(b)*(one-theta)*( omega*EVc_bu(ia) + (one-omega)*EVc_bu(ia+1) )
                        
                    	end if

                	end if

                end do

		    end do
	    end do
    	!$omp end parallel do

    end do
    
!!Solving for Decision Rules 
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!Set the Initial Distribution of Households

    !!!! Initialize Distribution of Households !!!!
    if ( iteration < 1 ) then
	    allocate(mu(ni,na,nj,ns))
	    allocate(mur(nb,na,nR))
    end if
    mu = zero; mur = zero

    !!!! Distribution of Households at Age 1 !!!!
    do i = 1,ni
	    e = e_index(i); n = n_index(i); b = b_index(i)
	    if ( n==1) then
		    mu(i,nd,1,1) = pi_erg(e)*phi_b(b)*phi*phi_e
	    end if
	    if ( n==2 ) then
	    	mu(i,nd,1,1) = pi_erg(e)*phi_b(b)*phi*(one-phi_e)
	    end if
    end do
    
    !!!! Future Distribution of Households - Working Ages !!!!
    iz = 1; im=25
    pi_i = pi_ig
    do j = 1,(nj-1)
        do a = 1,na
            do i = 1,ni
            
                !!!! Good Credit Status !!!!
                s = 1
                muval = mu(i,a,j,s)
                if ( muval > precision_mu ) then
                    gval = g(i,a,j,s,im)
    				dval = d(i,a,j,s,im)
	    			hval = h(i,a,j,s,im)
                    call interpolation(agrid, gval, 1 ,na, na ,ia, omega)
                
                    do i2 = 1,ni
                        if ( hval > (one-precision) ) then
                            if ( dval > (one-precision) ) then
                                mu(i2,ia,j+1,2) = mu(i2,ia,j+1,2) + omega*pi_i(i,i2,1)*muval
				    			mu(i2,ia+1,j+1,2) = mu(i2,ia+1,j+1,2) + (one-omega)*pi_i(i,i2,1)*muval
                            elseif ( dval < precision ) then
                                mu(i2,ia,j+1,1) = mu(i2,ia,j+1,1) + omega*pi_i(i,i2,1)*muval
							    mu(i2,ia+1,j+1,1) = mu(i2,ia+1,j+1,1) + (one-omega)*pi_i(i,i2,1)*muval
                            else
                                mu(i2,ia,j+1,1) = mu(i2,ia,j+1,1) + (one-dval)*omega*pi_i(i,i2,1)*muval
		    					mu(i2,ia+1,j+1,1) = mu(i2,ia+1,j+1,1) + (one-dval)*(one-omega)*pi_i(i,i2,1)*muval
			    				mu(i2,ia,j+1,2) = mu(i2,ia,j+1,2) + dval*omega*pi_i(i,i2,1)*muval
				    			mu(i2,ia+1,j+1,2) = mu(i2,ia+1,j+1,2) + dval*(one-omega)*pi_i(i,i2,1)*muval
                            end if
                        elseif ( hval < precision ) then
                            if ( dval > (one-precision) ) then
                                mu(i2,ia,j+1,2) = mu(i2,ia,j+1,2) + omega*pi_i(i,i2,2)*muval
			    				mu(i2,ia+1,j+1,2) = mu(i2,ia+1,j+1,2) + (one-omega)*pi_i(i,i2,2)*muval
                            elseif ( dval < precision ) then
                                mu(i2,ia,j+1,1) = mu(i2,ia,j+1,1) + omega*pi_i(i,i2,2)*muval
						    	mu(i2,ia+1,j+1,1) = mu(i2,ia+1,j+1,1) + (one-omega)*pi_i(i,i2,2)*muval
                            else
                                mu(i2,ia,j+1,1) = mu(i2,ia,j+1,1) + (one-dval)*omega*pi_i(i,i2,2)*muval
		    					mu(i2,ia+1,j+1,1) = mu(i2,ia+1,j+1,1) + (one-dval)*(one-omega)*pi_i(i,i2,2)*muval
			    				mu(i2,ia,j+1,2) = mu(i2,ia,j+1,2) + dval*omega*pi_i(i,i2,2)*muval
				    			mu(i2,ia+1,j+1,2) = mu(i2,ia+1,j+1,2) + dval*(one-omega)*pi_i(i,i2,2)*muval
                            end if
                        else
                            if ( dval > (one-precision) ) then
                                mu(i2,ia,j+1,2) = mu(i2,ia,j+1,2) + hval*omega*pi_i(i,i2,1)*muval
			    				mu(i2,ia+1,j+1,2) = mu(i2,ia+1,j+1,2) + hval*(one-omega)*pi_i(i,i2,1)*muval
				    			mu(i2,ia,j+1,2) = mu(i2,ia,j+1,2) + (one-hval)*omega*pi_i(i,i2,2)*muval
					    		mu(i2,ia+1,j+1,2) = mu(i2,ia+1,j+1,2) + (one-hval)*(one-omega)*pi_i(i,i2,2)*muval
                            elseif ( dval < precision ) then
                                mu(i2,ia,j+1,1) = mu(i2,ia,j+1,1) + hval*omega*pi_i(i,i2,1)*muval
    							mu(i2,ia+1,j+1,1) = mu(i2,ia+1,j+1,1) + hval*(one-omega)*pi_i(i,i2,1)*muval
	    						mu(i2,ia,j+1,1) = mu(i2,ia,j+1,1) + (one-hval)*omega*pi_i(i,i2,2)*muval
                                mu(i2,ia+1,j+1,1) = mu(i2,ia+1,j+1,1) + (one-hval)*(one-omega)*pi_i(i,i2,2)*muval
                            else
                                mu(i2,ia,j+1,1) = mu(i2,ia,j+1,1) + hval*(one-dval)*omega*pi_i(i,i2,1)*muval
					    		mu(i2,ia+1,j+1,1) = mu(i2,ia+1,j+1,1) + hval*(one-dval)*(one-omega)*pi_i(i,i2,1)*muval
						    	mu(i2,ia,j+1,1) = mu(i2,ia,j+1,1) + (one-hval)*(one-dval)*omega*pi_i(i,i2,2)*muval
							    mu(i2,ia+1,j+1,1) = mu(i2,ia+1,j+1,1) + (one-hval)*(one-dval)*(one-omega)*pi_i(i,i2,2)*muval
                    
    							mu(i2,ia,j+1,2) = mu(i2,ia,j+1,2) + hval*dval*omega*pi_i(i,i2,1)*muval
	    						mu(i2,ia+1,j+1,2) = mu(i2,ia+1,j+1,2) + hval*dval*(one-omega)*pi_i(i,i2,1)*muval
		    					mu(i2,ia,j+1,2) = mu(i2,ia,j+1,2) + (one-hval)*dval*omega*pi_i(i,i2,2)*muval
			    				mu(i2,ia+1,j+1,2) = mu(i2,ia+1,j+1,2) + (one-hval)*dval*(one-omega)*pi_i(i,i2,2)*muval
                            end if
                        end if
                    
                    end do
                end if
            
                !!!! Bad Credit Status !!!!
                s = 2
                muval = mu(i,a,j,s)
                if ( muval > precision_mu ) then
                    gval = g(i,a,j,s,im)
                    dval = d(i,a,j,s,im)
	    			hval = h(i,a,j,s,im)
                    call interpolation(agrid, gval, 1, na, na, ia, omega)
                
    				do i2 = 1,ni
                        if ( hval > (one-precision) ) then
		    				mu(i2,ia,j+1,1) = mu(i2,ia,j+1,1) + theta*omega*pi_i(i,i2,1)*muval
			    			mu(i2,ia+1,j+1,1) = mu(i2,ia+1,j+1,1) + theta*(one-omega)*pi_i(i,i2,1)*muval
                            mu(i2,ia,j+1,2) = mu(i2,ia,j+1,2) + (one-theta)*omega*pi_i(i,i2,1)*muval
					    	mu(i2,ia+1,j+1,2) = mu(i2,ia+1,j+1,2) + (one-theta)*(one-omega)*pi_i(i,i2,1)*muval
                        elseif ( hval < precision ) then
                            mu(i2,ia,j+1,1) = mu(i2,ia,j+1,1) + theta*omega*pi_i(i,i2,2)*muval
		    				mu(i2,ia+1,j+1,1) = mu(i2,ia+1,j+1,1) + theta*(one-omega)*pi_i(i,i2,2)*muval
                            mu(i2,ia,j+1,2) = mu(i2,ia,j+1,2) + (one-theta)*omega*pi_i(i,i2,2)*muval
				    		mu(i2,ia+1,j+1,2) = mu(i2,ia+1,j+1,2) + (one-theta)*(one-omega)*pi_i(i,i2,2)*muval
                        else
                            mu(i2,ia,j+1,1) = mu(i2,ia,j+1,1) + hval*theta*omega*pi_i(i,i2,1)*muval
    						mu(i2,ia+1,j+1,1) = mu(i2,ia+1,j+1,1) + hval*theta*(one-omega)*pi_i(i,i2,1)*muval
	    					mu(i2,ia,j+1,1) = mu(i2,ia,j+1,1) + (one-hval)*theta*omega*pi_i(i,i2,2)*muval
		    				mu(i2,ia+1,j+1,1) = mu(i2,ia+1,j+1,1) + (one-hval)*theta*(one-omega)*pi_i(i,i2,2)*muval
                    
    						mu(i2,ia,j+1,2) = mu(i2,ia,j+1,2) + hval*(one-theta)*omega*pi_i(i,i2,1)*muval
	    					mu(i2,ia+1,j+1,2) = mu(i2,ia+1,j+1,2) + hval*(one-theta)*(one-omega)*pi_i(i,i2,1)*muval
		    				mu(i2,ia,j+1,2) = mu(i2,ia,j+1,2) + (one-hval)*(one-theta)*omega*pi_i(i,i2,2)*muval
			    			mu(i2,ia+1,j+1,2) = mu(i2,ia+1,j+1,2) + (one-hval)*(one-theta)*(one-omega)*pi_i(i,i2,2)*muval
                        end if
                    end do
                end if

            end do
        end do
    end do
    
    !!!! Distribution of Households - Last Working Age !!!!
    j = nj
    do a = 1,na
        do i = 1,ni
            do s = 1,ns
			    muval = mu(i,a,j,s)
    			if ( muval > precision_mu ) then
                    b = b_index(i)
		    		gval = g(i,a,j,s,im)
                    call interpolation(agrid, gval, 1, na, na, ia, omega)
				    mur(b,ia,1) = mur(b,ia,1) + omega*muval
    				mur(b,ia+1,1) = mur(b,ia+1,1) + (one-omega)*muval
                end if
            
            end do
        end do
    end do

    !!!! Distribution of Households - Retirement Ages !!!!
    do j = 1,(nR-1)
        do a = 1,na
            do b = 1,nb
                muval = mur(b,a,j)
                if ( muval > precision_mu ) then
	    			gval = gr(b,a,j,im)
                    call interpolation(agrid, gval, 1, na, na, ia, omega)
                    mur(b,ia,j+1) = mur(b,ia,j+1) + omega*muval
                    mur(b,ia+1,j+1) = mur(b,ia+1,j+1) + (one-omega)*muval
                end if
            
            end do
        end do
    end do

    !!!! Save the Value Functions to Calculate Welfare - Age 1 !!!!
    if ( switch_results == 1 ) then
        do i = 1,ni
            e = e_index(i); n = n_index(i); b = b_index(i)
            open(unit=8, file='Welfare_Results.txt')
            write(unit=8, fmt='(I2, 3x, I1, 3x, I1, 3x, f16.13, 3x, f13.8, 3x, f13.8)') e, n, b, mu(i,nd,1,1), V(i,nd,1,1,im), V_c(i,nd,1,1,im)
        end do
    end if

!!Set the Initial Distribution of Households
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!Aggregate Simulation

    !!!! Allocate Large Arrays !!!!
    if ( iteration < 1 ) then
    	allocate(mu_next(ni,na,nj,ns))
	    allocate(mu_good(ni,na,nj,np))
	    allocate(mu_bad(ni,na,nj,np))
	    allocate(mur_next(nb,na,nR))
    end if

    !!!! Simulate Model Economy !!!!
    Dpop_P_mean = zero
    do t = 1,nt
	    iz = iz_sim(t)
        if ( iz == 1 ) then
            pi_i = pi_ig
        else
            pi_i = pi_ib
        endif

    	!!!! Calculate Aggregate Capital and Labor !!!!
	    aggK = zero; aggL = zero
        do i = 1,ni
            e = e_index(i); n = n_index(i)
            if ( n == 1 ) then
                do j = 1,nj
                    do a = 1,na
                        do s = 1,ns
                            muval = mu(i,a,j,s)
                            aggL = aggL + egrid(e)*gamma_j(j)*muval
                        end do
                    end do
                end do
            end if
        end do
        aggK = K_sim(t)
    
        !!!! Calculate Equilibrium Prices !!!!
        rate = zgrid(iz)*alpha*( aggK**(alpha-one) )*( aggL**(one-alpha) ) - delta
        wage = zgrid(iz)*(one-alpha)*( aggK**alpha )*( aggL**(-alpha) )
        call interpolation(kgrid, aggK, 1, nk, nk, ik, omega_k)
    	call interpolation(lgrid, aggL, 1, nl, nl, il, omega_l)
        do m = 1,nm
		    z = z_index(m); k = k_index(m); l = l_index(m)
    		if ( k == ik .and. l == il .and. z == iz ) then
	    		im1 = m
		    elseif ( k == ik+1 .and. l == il .and. z == iz ) then
    			im2 = m
	    	elseif ( k == ik .and. l == il+1 .and. z == iz ) then
		    	im3 = m
    		elseif ( k == ik+1 .and. l == il+1 .and. z == iz ) then
	    		im4 = m
		    end if
        end do
        rate_save = omega_k*omega_l*E_r(im1) + (one-omega_k)*omega_l*E_r(im2) + omega_k*(one-omega_l)*E_r(im3) + (one-omega_k)*(one-omega_l)*E_r(im4)
    
        !!!! Calculate Aggregate Future Capital !!!!
        aggKp = zero
        do j = 1,nj
            do i = 1,ni
                s = 1
                do a = 1,na
                    muval = mu(i,a,j,s)
                    gval = omega_k*omega_l*g(i,a,j,s,im1) + (one-omega_k)*omega_l*g(i,a,j,s,im2) + omega_k*(one-omega_l)*g(i,a,j,s,im3) + (one-omega_k)*(one-omega_l)*g(i,a,j,s,im4)
                    if ( gval < zero ) then
                        qval = omega_k*omega_l*q(i,a,j,im1) + (one-omega_k)*omega_l*q(i,a,j,im2) + omega_k*(one-omega_l)*q(i,a,j,im3) + (one-omega_k)*(one-omega_l)*q(i,a,j,im4)
                        aggKp = aggKp + qval*gval*muval
                    else
                        aggKp = aggKp + gval*muval/(one+rate_save)
                    end if
                end do
            
                s = 2
                do a = nd,na
                    muval = mu(i,a,j,s)
                    gval = omega_k*omega_l*g(i,a,j,s,im1) + (one-omega_k)*omega_l*g(i,a,j,s,im2) + omega_k*(one-omega_l)*g(i,a,j,s,im3) + (one-omega_k)*(one-omega_l)*g(i,a,j,s,im4)
                    aggKp = aggKp + gval*muval/(one+rate_save)
                end do
            end do
        
            if ( j <= nR ) then
                do b = 1,nb
                    do a = nd,na
                        muval = mur(b,a,j)
                        gval = omega_k*omega_l*gr(b,a,j,im1) + (1.0_rp-omega_k)*omega_l*gr(b,a,j,im2) + omega_k*(1.0_rp-omega_l)*gr(b,a,j,im3) + (1.0_rp-omega_k)*(1.0_rp-omega_l)*gr(b,a,j,im4)
                        aggKp = aggKp + gval*muval/(one+rate_save)
                    end do
                end do
            end if
        end do
    
        !!!! Collect Results !!!!
	    if ( switch_results == 1 ) then
        
            aggY = zero; aggC = zero; aggI = zero; aggG = zero; aggD = zero; aggB = zero; avg_earn = zero
		    Dpop = zero; Epop = zero; Upop = zero; Npop = zero; aggXi = zero
            Dpop_E = zero; Dpop_U = zero; Dpop_P = zero
            avg_q = zero; avg_limit = zero; pop_borrow = zero; limit_good = zero; pop_good = zero
            UI_exp = zero; ret_exp = zero; T_exp = zero; tax_rev = zero
            aggA = zero; agg_discharge = zero; profits=zero
        
            do a = 1,na
                do j = 1,nj
		    		do i = 1,ni
			    		e = e_index(i); n = n_index(i); b = b_index(i)
                        do s = 1,ns
                            muval = mu(i,a,j,s)
                            cval = omega_k*omega_l*c(i,a,j,s,im1) + (one-omega_k)*omega_l*c(i,a,j,s,im2) + omega_k*(one-omega_l)*c(i,a,j,s,im3) + (one-omega_k)*(one-omega_l)*c(i,a,j,s,im4)
	    					gval = omega_k*omega_l*g(i,a,j,s,im1) + (one-omega_k)*omega_l*g(i,a,j,s,im2) + omega_k*(one-omega_l)*g(i,a,j,s,im3) + (one-omega_k)*(one-omega_l)*g(i,a,j,s,im4)
		    				hval = omega_k*omega_l*h(i,a,j,s,im1) + (one-omega_k)*omega_l*h(i,a,j,s,im2) + omega_k*(one-omega_l)*h(i,a,j,s,im3) + (one-omega_k)*(one-omega_l)*h(i,a,j,s,im4)
			    			dval = omega_k*omega_l*d(i,a,j,s,im1) + (one-omega_k)*omega_l*d(i,a,j,s,im2) + omega_k*(one-omega_l)*d(i,a,j,s,im3) + (one-omega_k)*(one-omega_l)*d(i,a,j,s,im4)
                            blim = omega_k*omega_l*abar(i,a,j,im1) + (one-omega_k)*omega_l*abar(i,a,j,im2) + omega_k*(one-omega_l)*abar(i,a,j,im3) + (one-omega_k)*(one-omega_l)*abar(i,a,j,im4)
                        
                            if ( muval > precision_mu ) then
                            
                                !!!! Aggregate Quantities !!!!
                                aggC = aggC + cval*muval
                                aggB = aggB + dval*muval
                                T_exp = T_exp + y_bar*muval
                                if ( a < nd ) then
                                    aggD = aggD + agrid(a)*muval
                                    Dpop = Dpop + muval
                                    agg_discharge = agg_discharge + dval*agrid(a)*muval
                                    !blim = omega_k*omega_l*abar(i,a,j,im1) + (one-omega_k)*omega_l*abar(i,a,j,im2) + omega_k*(one-omega_l)*abar(i,a,j,im3) + (one-omega_k)*(one-omega_l)*abar(i,a,j,im4)
                                    avg_limit = avg_limit + blim*muval
                                else
                                    aggA = aggA + agrid(a)*muval
                                end if
                            
                                !!!! Results - Employment States !!!!
                                if ( n == 1 ) then
		    						Epop = Epop + muval
                                    aggXi = aggXi + (1.0_rp-hval)*muval
                                    aggXi = aggXi + xi(iz)*hval*muval
					    			avg_earn = avg_earn + wage*egrid(e)*gamma_j(j)*muval
						    		tax_rev = tax_rev + tau*wage*egrid(e)*gamma_j(j)*muval
							    	if ( a < nd ) then
									    Dpop_E = Dpop_E + muval
                                    end if
	    						else
		    						Upop = Upop + hval*muval
			    					Npop = Npop + (one-hval)*muval
				    				if ( a < nd ) then
	    								Dpop_U = Dpop_U + hval*muval
			    					end if
                                end if

                                !!!! UI Expenditures !!!!
	    						if ( n == 2 ) then 
		    						if ( upsilon_r(iz)*wage*egrid(e)*gamma_j(j) >= upsilon_bar(iz) ) then
			    						y_u = upsilon_bar(iz)
				    				else
					    				y_u = upsilon_r(iz)*wage*egrid(e)*gamma_j(j)
						    		end if
							    	UI_exp = UI_exp + y_u*muval
                                end if
                            
                                !!!! Borrowing and Saving !!!!
                                if ( gval < zero ) then
                                    qval = omega_k*omega_l*q(i,a,j,im1) + (one-omega_k)*omega_l*q(i,a,j,im2) + omega_k*(one-omega_l)*q(i,a,j,im3) + (one-omega_k)*(one-omega_l)*q(i,a,j,im4)
                                    pop_borrow = pop_borrow + muval
                                    avg_q = avg_q + qval*muval
                                end if
                                
                                !!!! Borrowing Limits !!!!
                                if ( s == 1 ) then
                                    limit_good = limit_good + blim*muval
                                    pop_good = pop_good + muval
                                end if

						    end if

					    end do					

                    end do

                    !!!! Results for Retired !!!!
	    			if ( j <= nR ) then
		    			do b = 1,nb

                            !!!! Aggregate Quantities !!!!
	    					cval = omega_k*omega_l*cr(b,a,j,im1) + (1.0_rp-omega_k)*omega_l*cr(b,a,j,im2) + omega_k*(1.0_rp-omega_l)*cr(b,a,j,im3) + (1.0_rp-omega_k)*(1.0_rp-omega_l)*cr(b,a,j,im4)
		    				gval = omega_k*omega_l*gr(b,a,j,im1) + (1.0_rp-omega_k)*omega_l*gr(b,a,j,im2) + omega_k*(1.0_rp-omega_l)*gr(b,a,j,im3) + (1.0_rp-omega_k)*(1.0_rp-omega_l)*gr(b,a,j,im4)
			    			aggC = aggC + cval*mur(b,a,j)

    						ret_exp = ret_exp + y_r*mur(b,a,j)
                            if ( a > nd ) then
                        		aggA = aggA + agrid(a)*mur(b,a,j)
                    		end if

    					end do
	    			end if

    			end do

            end do
            aggXi = aggXi/Epop
            avg_earn = avg_earn/Epop
            avg_limit = avg_limit/Dpop
            limit_good = limit_good/pop_good
            Dpop = Dpop/sum(mu)
            Dpop_P = (Dpop_E+Dpop_U)/(Epop+Upop)
            if ( t > 300 ) then
                Dpop_P_mean = Dpop_P_mean + Dpop_P
            end if
        
		    Npop = Npop + sum(mur)
            aggI = aggKp - (one-delta)*aggK
		    aggY = zgrid(iz)*((aggK)**(alpha))*(aggL**(one-alpha))
            profits = (one+rate)*aggK - aggA - aggD + agg_discharge + iota(iz)*aggD
            aggG = tax_rev - UI_exp - ret_exp - T_exp + profits
            avg_q = avg_q/pop_borrow

	    end if

	    !!!! Calculate Distribution of Households !!!!
	    ip = 1
	    mu_bad = zero; mu_good = zero; mur_next = zero; mu_next = zero
	    mu_next(:,:,1,:) = mu(:,:,1,:)

	    !!!! Future Distribution of Households - Working Ages !!!!
        !$omp parallel do private(j, a, i, b, s, i2, ip, ia, omega, gval, hval, dval, muval)
	    do j = 1,(nj-1)
		    do a = 1,na
			    do i = 1,ni
				    !$ ip = omp_get_thread_num()
                    !$ ip = ip + 1
                
				    !!!! Good Credit Status !!!!
				    s = 1
				    muval = mu(i,a,j,s)
				    if ( muval > precision_mu ) then
					    gval = omega_k*omega_l*g(i,a,j,s,im1) + (one-omega_k)*omega_l*g(i,a,j,s,im2) + omega_k*(one-omega_l)*g(i,a,j,s,im3) + (one-omega_k)*(one-omega_l)*g(i,a,j,s,im4)
					    hval = omega_k*omega_l*h(i,a,j,s,im1) + (one-omega_k)*omega_l*h(i,a,j,s,im2) + omega_k*(one-omega_l)*h(i,a,j,s,im3) + (one-omega_k)*(one-omega_l)*h(i,a,j,s,im4)
					    dval = omega_k*omega_l*d(i,a,j,s,im1) + (one-omega_k)*omega_l*d(i,a,j,s,im2) + omega_k*(one-omega_l)*d(i,a,j,s,im3) + (one-omega_k)*(one-omega_l)*d(i,a,j,s,im4)
                        call interpolation(agrid, gval, 1 ,na, na ,ia, omega)
					    do i2 = 1,ni
						    if ( hval > (one-precision) ) then
							    if ( dval > (one-precision) ) then
								    mu_bad(i2,nd,j+1,ip) = mu_bad(i2,nd,j+1,ip) + omega*pi_i(i,i2,1)*muval
    							elseif ( dval < precision ) then
	    							mu_good(i2,ia,j+1,ip) = mu_good(i2,ia,j+1,ip) + omega*pi_i(i,i2,1)*muval
		    						mu_good(i2,ia+1,j+1,ip) = mu_good(i2,ia+1,j+1,ip) + (one-omega)*pi_i(i,i2,1)*muval
			    				else
				    				mu_good(i2,ia,j+1,ip) = mu_good(i2,ia,j+1,ip) + (one-dval)*omega*pi_i(i,i2,1)*muval
					    			mu_good(i2,ia+1,j+1,ip) = mu_good(i2,ia+1,j+1,ip) + (one-dval)*(one-omega)*pi_i(i,i2,1)*muval
						    		mu_bad(i2,nd,j+1,ip) = mu_bad(i2,nd,j+1,ip) + dval*omega*pi_i(i,i2,1)*muval
							    end if
    						elseif ( hval < precision ) then
	    						if ( dval > (one-precision) ) then
		    						mu_bad(i2,nd,j+1,ip) = mu_bad(i2,nd,j+1,ip) + omega*pi_i(i,i2,2)*muval
				    			elseif ( dval < precision ) then
					    			mu_good(i2,ia,j+1,ip) = mu_good(i2,ia,j+1,ip) + omega*pi_i(i,i2,2)*muval
						    		mu_good(i2,ia+1,j+1,ip) = mu_good(i2,ia+1,j+1,ip) + (one-omega)*pi_i(i,i2,2)*muval
							    else
    								mu_good(i2,ia,j+1,ip) = mu_good(i2,ia,j+1,ip) + (one-dval)*omega*pi_i(i,i2,2)*muval
	    							mu_good(i2,ia+1,j+1,ip) = mu_good(i2,ia+1,j+1,ip) + (one-dval)*(one-omega)*pi_i(i,i2,2)*muval
		    						mu_bad(i2,nd,j+1,ip) = mu_bad(i2,nd,j+1,ip) + dval*omega*pi_i(i,i2,2)*muval
				    			end if
					    	else
    							if ( dval > (one-precision) ) then
	    							mu_bad(i2,nd,j+1,ip) = mu_bad(i2,nd,j+1,ip) + hval*omega*pi_i(i,i2,1)*muval
			    					mu_bad(i2,nd,j+1,ip) = mu_bad(i2,nd,j+1,ip) + (one-hval)*omega*pi_i(i,i2,2)*muval
					    		elseif ( dval < precision ) then
						    		mu_good(i2,ia,j+1,ip) = mu_good(i2,ia,j+1,ip) + hval*omega*pi_i(i,i2,1)*muval
							    	mu_good(i2,ia+1,j+1,ip) = mu_good(i2,ia+1,j+1,ip) + hval*(one-omega)*pi_i(i,i2,1)*muval
								    mu_good(i2,ia,j+1,ip) = mu_good(i2,ia,j+1,ip) + (one-hval)*omega*pi_i(i,i2,2)*muval
    								mu_good(i2,ia+1,j+1,ip) = mu_good(i2,ia+1,j+1,ip) + (one-hval)*(one-omega)*pi_i(i,i2,2)*muval
	    						else
								    mu_good(i2,ia,j+1,ip) = mu_good(i2,ia,j+1,ip) + hval*(one-dval)*omega*pi_i(i,i2,1)*muval
    								mu_good(i2,ia+1,j+1,ip) = mu_good(i2,ia+1,j+1,ip) + hval*(one-dval)*(one-omega)*pi_i(i,i2,1)*muval
	    							mu_bad(i2,nd,j+1,ip) = mu_bad(i2,nd,j+1,ip) + hval*dval*omega*pi_i(i,i2,1)*muval
                            
    								mu_good(i2,ia,j+1,ip) = mu_good(i2,ia,j+1,ip) + (one-hval)*(one-dval)*omega*pi_i(i,i2,2)*muval
	    							mu_good(i2,ia+1,j+1,ip) = mu_good(i2,ia+1,j+1,ip) + (one-hval)*(one-dval)*(one-omega)*pi_i(i,i2,2)*muval
		    						mu_bad(i2,nd,j+1,ip) = mu_bad(i2,nd,j+1,ip) + (one-hval)*dval*omega*pi_i(i,i2,2)*muval
				    			end if
					    	end if
                    
					    end do
				    end if
            
    				!!!! Bad Credit Status !!!!
	    			s = 2
		    		muval = mu(i,a,j,s)
			    	if ( muval > precision_mu ) then
    					gval = omega_k*omega_l*g(i,a,j,s,im1) + (one-omega_k)*omega_l*g(i,a,j,s,im2) + omega_k*(one-omega_l)*g(i,a,j,s,im3) + (one-omega_k)*(one-omega_l)*g(i,a,j,s,im4)
	    				hval = omega_k*omega_l*h(i,a,j,s,im1) + (one-omega_k)*omega_l*h(i,a,j,s,im2) + omega_k*(one-omega_l)*h(i,a,j,s,im3) + (one-omega_k)*(one-omega_l)*h(i,a,j,s,im4)
		    			dval = omega_k*omega_l*d(i,a,j,s,im1) + (one-omega_k)*omega_l*d(i,a,j,s,im2) + omega_k*(one-omega_l)*d(i,a,j,s,im3) + (one-omega_k)*(one-omega_l)*d(i,a,j,s,im4)
                        call interpolation(agrid, gval, 1, na, na, ia, omega)
				    	do i2 = 1,ni
    						if ( hval > (one-precision) ) then
	    						mu_good(i2,ia,j+1,ip) = mu_good(i2,ia,j+1,ip) + theta*omega*pi_i(i,i2,1)*muval
		    					mu_good(i2,ia+1,j+1,ip) = mu_good(i2,ia+1,j+1,ip) + theta*(one-omega)*pi_i(i,i2,1)*muval
			    				mu_bad(i2,ia,j+1,ip) = mu_bad(i2,ia,j+1,ip) + (one-theta)*omega*pi_i(i,i2,1)*muval
				    			mu_bad(i2,ia+1,j+1,ip) = mu_bad(i2,ia+1,j+1,ip) + (one-theta)*(one-omega)*pi_i(i,i2,1)*muval
					    	elseif ( hval < precision ) then
    							mu_good(i2,ia,j+1,ip) = mu_good(i2,ia,j+1,ip) + theta*omega*pi_i(i,i2,2)*muval
	    						mu_good(i2,ia+1,j+1,ip) = mu_good(i2,ia+1,j+1,ip) + theta*(one-omega)*pi_i(i,i2,2)*muval
		    					mu_bad(i2,ia,j+1,ip) = mu_bad(i2,ia,j+1,ip) + (one-theta)*omega*pi_i(i,i2,2)*muval
			    				mu_bad(i2,ia+1,j+1,ip) = mu_bad(i2,ia+1,j+1,ip) + (one-theta)*(one-omega)*pi_i(i,i2,2)*muval
				    		else
					    		mu_good(i2,ia,j+1,ip) = mu_good(i2,ia,j+1,ip) + hval*theta*omega*pi_i(i,i2,1)*muval
						    	mu_good(i2,ia+1,j+1,ip) = mu_good(i2,ia+1,j+1,ip) + hval*theta*(one-omega)*pi_i(i,i2,1)*muval
							    mu_good(i2,ia,j+1,ip) = mu_good(i2,ia,j+1,ip) + (one-hval)*theta*omega*pi_i(i,i2,2)*muval
    							mu_good(i2,ia+1,j+1,ip) = mu_good(i2,ia+1,j+1,ip) + (one-hval)*theta*(one-omega)*pi_i(i,i2,2)*muval
                    
    							mu_bad(i2,ia,j+1,ip) = mu_bad(i2,ia,j+1,ip) + hval*(one-theta)*omega*pi_i(i,i2,1)*muval
	    						mu_bad(i2,ia+1,j+1,ip) = mu_bad(i2,ia+1,j+1,ip) + hval*(one-theta)*(one-omega)*pi_i(i,i2,1)*muval
		    					mu_bad(i2,ia,j+1,ip) = mu_bad(i2,ia,j+1,ip) + (one-hval)*(one-theta)*omega*pi_i(i,i2,2)*muval
			    				mu_bad(i2,ia+1,j+1,ip) = mu_bad(i2,ia+1,j+1,ip) + (one-hval)*(one-theta)*(one-omega)*pi_i(i,i2,2)*muval
				    		end if
					    end do
				    end if

			    end do
		    end do
        end do
        !$omp end parallel do
    
    	!!!! Distribution of Households - Last Working Age !!!!
	    j = nj
    	do a = 1,na
	    	do i = 1,ni
		    	do s = 1,ns     
			    	muval = mu(i,a,j,s)
				    if ( muval > precision_mu ) then
    					gval = omega_k*omega_l*g(i,a,j,s,im1) + (one-omega_k)*omega_l*g(i,a,j,s,im2) + omega_k*(one-omega_l)*g(i,a,j,s,im3) + (one-omega_k)*(one-omega_l)*g(i,a,j,s,im4)
	    				b = b_index(i)
                        call interpolation(agrid, gval, 1 ,na, na ,ia, omega)
			    		mur_next(b,ia,1) = mur_next(b,ia,1) + omega*muval
				    	mur_next(b,ia+1,1) = mur_next(b,ia+1,1) + (one-omega)*muval
				    end if
			    end do
		    end do
	    end do

    	!!!! Distribution of Households - Retirement Ages !!!!
	    do j = 1,(nR-1)
		    do a = 1,na
			    do b = 1,nb
				    muval = mur(b,a,j)
    				if ( muval > precision_mu ) then
	    				gval = omega_k*omega_l*gr(b,a,j,im1) + (one-omega_k)*omega_l*gr(b,a,j,im2) + omega_k*(one-omega_l)*gr(b,a,j,im3) + (one-omega_k)*(one-omega_l)*gr(b,a,j,im4)
                        call interpolation(agrid, gval, 1, na, na, ia, omega)
			    		mur_next(b,ia,j+1) = mur_next(b,ia,j+1) + omega*muval
				    	mur_next(b,ia+1,j+1) = mur_next(b,ia+1,j+1) + (one-omega)*muval
				    end if
			    end do
		    end do
        end do

	    !!!! Solve for Future Distribution !!!!
        do j = 2,nj
            do a = 1,na
                do i = 1,ni
                    mu(i,a,j,1) = sum(mu_good(i,a,j,:))
                    mu(i,a,j,2) = sum(mu_bad(i,a,j,:))
                end do
            end do
        end do
        mur(:,:,:) = mur_next(:,:,:)

    	!!!! Save Simulation Results !!!!
        if ( switch_results == 1 ) then
            Z_sim(t) = zgrid(iz)
            lambda_sim(t) = lambda(iz)
            xi_sim(t) = aggXi
    		L_sim(t) = aggL
            K_sim(t+1) = aggKp
            Y_sim(t) = aggY
            C_sim(t) = aggC
		    I_sim(t) = aggI
            G_sim(t) = aggG
	    	D_sim(t) = aggD
		    B_sim(t) = aggB
    		U_sim(t) = 100.0_rp*Upop/(Epop+Upop)
	    	P_sim(t) = 100.0_rp*(Epop+Upop)/(Epop+Upop+Npop)
            r_sim(t) = rate
	    	w_sim(t) = wage
            M_sim(t) = 100.0_rp*( ( (2.0_rp-avg_q)**4.0_rp - 1.0_rp ) - ( (1.0_rp+rate_save)**4.0_rp - 1.0_rp ) )
            abar_sim(t) = limit_good
            ebar_sim(t) = avg_earn
        else
            Z_sim(t) = zgrid(iz)
            L_sim(t) = aggL
            K_sim(t+1) = aggKp
            r_sim(t) = rate
            w_sim(t) = wage
        end if

	    !!!! Save/Display Simulation Results !!!!
        if ( switch_results == 1 ) then
	    	print *, ""
		    print *, "-----------------------------------------"
            print *, ""
	    	print '(1x, a6, 3x, a9)', "Period", "Agg State"
		    print '(1x, I4, 6x, I4)', t, iz
            print *, ""
        
            open(unit=1, file='Exogenous_Fluctuations.txt')
            write(unit=1, fmt='(f10.6, f10.6, f10.6)') zgrid(iz), aggXi, lambda(iz)
        
            open(unit=2, file='Forecasting_Rules.txt')
            write(unit=2, fmt='(I1, 4x, f10.6, f10.6)') iz, aggK, aggL
        
            open(unit=3, file='Equilibrium_Prices.txt')
            write(unit=3, fmt='(f10.6, f10.6)') rate, wage
        
            open(unit=4, file='Simulation_Results1.txt')
            write(unit=4, fmt='(f10.6, f10.6, f10.6)') aggY, aggC, aggI
        
            open(unit=5, file='Simulation_Results2.txt')
            write(unit=5, fmt='(f10.6, f10.6, f11.7)') aggD, aggB, 100.0_rp*( ( (2.0_rp-avg_q)**4.0_rp - 1.0_rp ) - ( (1.0_rp+rate_save)**4.0_rp - 1.0_rp ) )
        
            open(unit=6, file='Simulation_Results3.txt')
            write(unit=6, fmt='(f10.6, f10.6, f10.6)') Upop/(Epop+Upop), (Epop+Upop)/(Epop+Upop+Npop), Dpop
        
            open(unit=7, file='Simulation_Results4.txt')
            write(unit=7, fmt='(f10.6, f10.6, f10.6)') aggG, avg_limit, limit_good
            
        end if

    end do

!! Aggregate Simulation
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!! End Outer Loop
    
    !!!! Length of Expansions and Recessions !!!!
    count_g = 0; count_b = 0
    do t = 301,(nt-1)
        iz = iz_sim(t)
        if ( iz == 2 ) then
            count_b = count_b + 1
        else
            count_g = count_g + 1
        end if
    end do

    !!!! Capital and Labor Moments !!!!
    avg_k = sum(K_sim(301:nt))/dble(nt-300)
    avg_l = sum(L_sim(301:nt))/dble(nt-300)
    sigma_K = zero
    sigma_L = zero
    do t = 301,nt
        sigma_K = sigma_K + (K_sim(t)-avg_K)**2.0_rp
        sigma_L = sigma_L + (L_sim(t)-avg_L)**2.0_rp
    end do
    sigma_K = sqrt( sigma_K/dble(nt-301-1) )
    sigma_L = sqrt( sigma_L/dble(nt-301-1) )
    
    !!!! Create Data for Regressions !!!!
    if ( iteration < 1 ) then
	    allocate( K_sim_g(count_g), L_sim_g(count_g), K_sim_b(count_b), L_sim_b(count_b), X_g(count_g,3), X_b(count_b,3) )
    end if

    count_g = 0; count_b = 0
    do t = 301,(nt-1)
        iz = iz_sim(t)
        if ( iz == 2 ) then
            count_b = count_b + 1
            X_b(count_b,1) = one
            X_b(count_b,2) = log( K_sim(t) )
            X_b(count_b,3) = log( L_sim(t) )
            K_sim_b(count_b) = log( K_sim(t+1) )
            L_sim_b(count_b) = log( L_sim(t+1) )
        else
            count_g = count_g + 1
            X_g(count_g,1) = one
            X_g(count_g,2) = log( K_sim(t) )
            X_g(count_g,3) = log( L_sim(t) )
            K_sim_g(count_g) = log( K_sim(t+1) )
            L_sim_g(count_g) = log( L_sim(t+1) )
        end if
    end do

    !!!! Calculate Regression Coefficients !!!!
    call regression_coefficients(count_g, 3, K_sim_g, X_g, b_g)
    nu_p_k0(1) = b_g(1)
    nu_p_k1(1) = b_g(2)
    nu_p_k2(1) = b_g(3)
    call regression_coefficients(count_g, 3, L_sim_g, X_g, b_g)
    nu_p_l0(1) = b_g(1)
    nu_p_l1(1) = b_g(3)
    nu_p_l2(1) = b_g(2)
    call regression_coefficients(count_b, 3, K_sim_b, X_b, b_b)
    nu_p_k0(2) = b_b(1)
    nu_p_k1(2) = b_b(2)
    nu_p_k2(2) = b_b(3)
    call regression_coefficients(count_b, 3, L_sim_b, X_b, b_b)
    nu_p_l0(2) = b_b(1)
    nu_p_l1(2) = b_b(3)
    nu_p_l2(2) = b_b(2)

    !!!! Save Forecasting Rules !!!!
    open(unit=22, file='Quantitative Results.txt')
    write(unit=22, fmt=*) "Scenario= ", switch_scenario
    write(unit=22, fmt=*) "Agg Vars= ", zgrid(2), xi(2), lambda(2)
    write(unit=22, fmt=*), "UI Vars= ", upsilon_d(2), upsilon_r(2), upsilon_bar(2)
    write(unit=22, fmt=*), ""
    write(unit=22, fmt=*), ""
    write(unit=22, fmt=*), "----------------------------------------------------------------------"
    write(unit=22, fmt=*), "---------------- Regression Coefficients - Original ------------------"
    write(unit=22, fmt=*), ""
    write(unit=22, fmt='(1x, a17, 3x, f12.8, 3x, f12.8 )'), "Capital Moments: ", avg_K, sigma_K
    write(unit=22, fmt='(1x, a15, 5x, f12.8, 3x, f12.8, 3x, f12.8 )'), "Capital (z_g): ", nu_k0(1), nu_k1(1), nu_k2(1)
    write(unit=22, fmt='(1x, a15, 5x, f12.8, 3x, f12.8, 3x, f12.8 )'), "Capital (z_b): ", nu_k0(2), nu_k1(2), nu_k2(2)
    write(unit=22, fmt=*), ""
    write(unit=22, fmt='(1x, a15, 5x, f12.8, 3x, f12.8 )'), "Labor Moments: ", avg_L, sigma_L
    write(unit=22, fmt='(1x, a13, 7x, f12.8, 3x, f12.8, 3x, f12.8 )'), "Labor (z_g): ", nu_l0(1), nu_l1(1), nu_l2(1)
    write(unit=22, fmt='(1x, a13, 7x, f12.8, 3x, f12.8, 3x, f12.8 )'), "Labor (z_b): ", nu_l0(2), nu_l1(2), nu_l2(2)
    write(unit=22, fmt=*), ""

    !!!! End Outer Loop !!!!
    if ( switch_scenario == 11 ) then
        precision_outer = 0.002_rp
    else
        precision_outer = 0.001_rp
    end if
    error = maxval( (/ abs(nu_k0(1)-nu_p_k0(1)), abs(nu_k1(1)-nu_p_k1(1)), abs(nu_k2(1)-nu_p_k2(1)), abs(nu_k0(2)-nu_p_k0(2)), abs(nu_k1(2)-nu_p_k1(2)), abs(nu_k2(2)-nu_p_k2(2)) /) )
    error = maxval( (/ error, abs(nu_l0(1)-nu_p_l0(1)), abs(nu_l1(1)-nu_p_l1(1)), abs(nu_l2(1)-nu_p_l2(1)), abs(nu_l0(2)-nu_p_l0(2)), abs(nu_l1(2)-nu_p_l1(2)), abs(nu_l2(2)-nu_p_l2(2)) /) )
    if ( iteration < max_iter ) then
        if ( error <= precision_outer ) then
            iteration = max_iter
            switch_results = 1
        else 
            iteration = iteration + 1
            if ( iteration < 10 ) then
                nu_k0(1) = 0.5_rp*nu_k0(1)+0.5_rp*nu_p_k0(1); nu_k1(1) = 0.5_rp*nu_k1(1)+0.5_rp*nu_p_k1(1); nu_k2(1) = 0.5_rp*nu_k2(1)+0.5_rp*nu_p_k2(1)
                nu_k0(2) = 0.5_rp*nu_k0(2)+0.5_rp*nu_p_k0(2); nu_k1(2) = 0.5_rp*nu_k1(2)+0.5_rp*nu_p_k1(2); nu_k2(2) = 0.5_rp*nu_k2(2)+0.5_rp*nu_p_k2(2)
                nu_l0(1) = 0.5_rp*nu_l0(1)+0.5_rp*nu_p_l0(1); nu_l1(1) = 0.5_rp*nu_l1(1)+0.5_rp*nu_p_l1(1); nu_l2(1) = 0.5_rp*nu_l2(1)+0.5_rp*nu_p_l2(1)
                nu_l0(2) = 0.5_rp*nu_l0(2)+0.5_rp*nu_p_l0(2); nu_l1(2) = 0.5_rp*nu_l1(2)+0.5_rp*nu_p_l1(2); nu_l2(2) = 0.5_rp*nu_l2(2)+0.5_rp*nu_p_l2(2)
            else
                nu_k0(1) = 0.75_rp*nu_k0(1)+0.25_rp*nu_p_k0(1); nu_k1(1) = 0.75_rp*nu_k1(1)+0.25_rp*nu_p_k1(1); nu_k2(1) = 0.75_rp*nu_k2(1)+0.25_rp*nu_p_k2(1)
                nu_k0(2) = 0.75_rp*nu_k0(2)+0.25_rp*nu_p_k0(2); nu_k1(2) = 0.75_rp*nu_k1(2)+0.25_rp*nu_p_k1(2); nu_k2(2) = 0.75_rp*nu_k2(2)+0.25_rp*nu_p_k2(2)
                nu_l0(1) = 0.75_rp*nu_l0(1)+0.25_rp*nu_p_l0(1); nu_l1(1) = 0.75_rp*nu_l1(1)+0.25_rp*nu_p_l1(1); nu_l2(1) = 0.75_rp*nu_l2(1)+0.25_rp*nu_p_l2(1)
                nu_l0(2) = 0.75_rp*nu_l0(2)+0.25_rp*nu_p_l0(2); nu_l1(2) = 0.75_rp*nu_l1(2)+0.25_rp*nu_p_l1(2); nu_l2(2) = 0.75_rp*nu_l2(2)+0.25_rp*nu_p_l2(2)
            end if
        end if
    else
        iteration = iteration + 1
    end if

    print *, "Outer Loop: ", iteration, error
    write(unit=22, fmt=*), "Outer Loop: ", iteration, error
    
end do
    
!! End Outer Loop
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!! Individual Simulation
    
!!!! Allocate Arrays for Individual Simulation !!!!
allocate(rand_epsilon(nj,nh,nsim), rand_unemp(nj,nh,nsim), rand_beta(nj,nh,nsim), rand_credit(nj,nh,nsim), rand_upsilon(nj,nh,nsim))
allocate(c_isim(nj,nh,nsim), h_isim(nj,nh,nsim), g_isim(nj,nh,nsim), d_isim(nj,nh,nsim), q_isim(nj,nh,nsim), inc_isim(nj,nh,nsim), abar_isim(nj,nh,nsim))
    
print *, ""
print *, ""
print *, "--------------------------------"
print *, " Simulate Individual Households "
print *, "--------------------------------"
print *, ""
    
!!!! Upload Aggregate State Space !!!!
open( unit=1, file="Forecasting_Rules.txt", status='old' )
do t = 1,nt
    read(1,*) z_sim(t), K_sim(t), L_sim(t)
end do

!!!! Random Productivity Shocks !!!!
call RANDOM_SEED( size = seedsize )
seed = (/ 516641, 2112222 /)
call RANDOM_SEED(put=seed)
call RANDOM_NUMBER(rand_epsilon)

!!!! Random Unemployment Shocks !!!!
call RANDOM_SEED(size = seedsize)
seed = (/ 616641, 3112222 /)
call RANDOM_SEED(put=seed)
call RANDOM_NUMBER(rand_unemp)
    
!!!! Random UI Shocks !!!!
call RANDOM_SEED(size = seedsize)
seed = (/ 716641, 4112222 /)
call RANDOM_SEED(put=seed)
call RANDOM_NUMBER(rand_upsilon)
    
!!!! Random Credit Shocks !!!!
call RANDOM_SEED(size = seedsize)
seed = (/ 816641, 4112222 /)
call RANDOM_SEED(put=seed)
call RANDOM_NUMBER(rand_credit)
    
!!!! Random Discount Factor Shocks !!!!
call RANDOM_SEED(size=seedsize)
seed = (/ 816641, 5112222 /)
call RANDOM_SEED(put=seed)
call RANDOM_NUMBER(rand_beta)
    
!!!! Initial Individual State Space !!!!
do t = 1,nsim
    do ih = 1,nh
        ia_sim(1,ih,t) = nd
        is_sim(1,ih,t) = 1
            
        call future_index(cum_pi_erg, rand_epsilon(1,ih,t), ne, ie_sim(1,ih,t))
            
        if ( rand_beta(1,ih,t) <= pi_beta ) then
            ib_sim(1,ih,t) = 1
        else
            ib_sim(1,ih,t) = 2
        end if
            
        if ( rand_unemp(1,ih,t) <= phi_e ) then
            in_sim(1,ih,t) = 1
        else
            in_sim(1,ih,t) = 2
        end if
    end do
end do
    
!!!! Begin Simulation !!!!
it = 0
mpc_mean = zero; count_mpc = zero
do t = 980,1188
    it = it+1
        
    !!!! Discretize Aggregate State !!!!
    iz = iz_sim(t); aggK = K_sim(t); aggL = L_sim(t)
    call interpolation(kgrid, aggK, 1, nk, nk, ik, omega_k)
    call interpolation(lgrid, aggL, 1, nl, nl, il, omega_l)
    rate = zgrid(iz)*alpha*(aggK**(alpha-one))*(aggL**(one-alpha)) - delta
    wage = zgrid(iz)*(one-alpha)*(aggK**alpha)*(aggL**(-alpha))
    do m = 1,nm
        z = z_index(m); k = k_index(m); l = l_index(m)
        if ( k == ik .and. l == il .and. z == iz ) then
            im1 = m
        elseif ( k == ik+1 .and. l == il .and. z == iz ) then
    	    im2 = m
	    elseif ( k == ik .and. l == il+1 .and. z == iz ) then
            im3 = m
        elseif ( k == ik+1 .and. l == il+1 .and. z == iz ) then
    	    im4 = m
	    end if
    end do
    rate_save = omega_k*omega_l*E_r(im1) + (one-omega_k)*omega_l*E_r(im2) + omega_k*(one-omega_l)*E_r(im3) + (one-omega_k)*(one-omega_l)*E_r(im4)
        
    do j = 1,nj
        do ih = 1,nh
            e = ie_sim(j,ih,it); n = in_sim(j,ih,it); b = ib_sim(j,ih,it); a = ia_sim(j,ih,it); s = is_sim(j,ih,it)
            i = i_index(e,n,b)
                
            if ( j <= it ) then
                
                !!!! Individual Decision Rules !!!!
                cval = omega_k*omega_l*c(i,a,j,s,im1) + (one-omega_k)*omega_l*c(i,a,j,s,im2) + omega_k*(one-omega_l)*c(i,a,j,s,im3) + (one-omega_k)*(one-omega_l)*c(i,a,j,s,im4)
                gval = omega_k*omega_l*g(i,a,j,s,im1) + (one-omega_k)*omega_l*g(i,a,j,s,im2) + omega_k*(one-omega_l)*g(i,a,j,s,im3) + (one-omega_k)*(one-omega_l)*g(i,a,j,s,im4)
                hval = omega_k*omega_l*h(i,a,j,s,im1) + (one-omega_k)*omega_l*h(i,a,j,s,im2) + omega_k*(one-omega_l)*h(i,a,j,s,im3) + (one-omega_k)*(one-omega_l)*h(i,a,j,s,im4)
                dval = omega_k*omega_l*d(i,a,j,s,im1) + (one-omega_k)*omega_l*d(i,a,j,s,im2) + omega_k*(one-omega_l)*d(i,a,j,s,im3) + (one-omega_k)*(one-omega_l)*d(i,a,j,s,im4)
                if ( s == 1 .and. gval < zero ) then
                    qval = omega_k*omega_l*q(i,a,j,im1) + (one-omega_k)*omega_l*q(i,a,j,im2) + omega_k*(one-omega_l)*q(i,a,j,im3) + (one-omega_k)*(one-omega_l)*q(i,a,j,im4)
                else
                qval = one/(one+rate_save)
                end if
                c_isim(j,ih,it) = cval
                g_isim(j,ih,it) = gval
                h_isim(j,ih,it) = hval
                d_isim(j,ih,it) = dval
                q_isim(j,ih,it) = qval
                if ( s == 1 ) then
                    abar_isim(j,ih,it) = omega_k*omega_l*abar(i,a,j,im1) + (one-omega_k)*omega_l*abar(i,a,j,im2) + omega_k*(one-omega_l)*abar(i,a,j,im3) + (one-omega_k)*(one-omega_l)*abar(i,a,j,im4)
                else
                    abar_isim(j,ih,it) = zero
                end if
                    
                !!!! Calculate MPC !!!!
                if ( t == 1170 ) then
                    a_dev = agrid(a) + 0.04_rp
                    call interpolation(agrid, a_dev, 1, na, na, ia ,omega_a)
                    c_dev = omega_a*( omega_k*omega_l*c(i,ia,j,s,im1) + (one-omega_k)*omega_l*c(i,ia,j,s,im2) + omega_k*(one-omega_l)*c(i,ia,j,s,im3) + (one-omega_k)*(one-omega_l)*c(i,ia,j,s,im4) ) + (one-omega_a)*( omega_k*omega_l*c(i,ia+1,j,s,im1) + (one-omega_k)*omega_l*c(i,ia+1,j,s,im2) + omega_k*(one-omega_l)*c(i,ia+1,j,s,im3) + (one-omega_k)*(one-omega_l)*c(i,ia+1,j,s,im4) )
                    count_mpc = count_mpc + one
                    mpc_mean = mpc_mean + (c_dev - cval)/0.04_rp
                end if
                    
                if ( j < nj ) then  
                    
                    !!!! Future Asset State !!!!
                    call interpolation(agrid, gval, 1, na, na, ia, omega_a) 
                    if ( omega_a <= 0.5_rp ) then
                        ia_sim(j+1,ih,it+1) = ia+1
                    else
                        ia_sim(j+1,ih,it+1) = ia
                    end if
                    
                    !!!! Future Employment State and Future Productivity State !!!!
                    if ( hval <= 0.5_rp ) then
                        in_sim(j+1,ih,it+1) = 3
                        ie_sim(j+1,ih,it+1) = e
                    else
                        if ( n == 1 ) then
                            if ( rand_unemp(j,ih,it) <= xi(iz) ) then
                                in_sim(j+1,ih,it+1) = 2
                                ie_sim(j+1,ih,it+1) = e
                            else
                                in_sim(j+1,ih,it+1) = 1
                                call future_index(cum_pi_e(e,:), rand_epsilon(j+1,ih,it+1), ne, ie_sim(j+1,ih,it+1))
                            end if
                        elseif (n == 2 ) then
                            if ( rand_unemp(j,ih,it) <= lambda(iz) ) then
                                in_sim(j+1,ih,it+1) = 1
                                call future_index(cum_pi_e(e,:), rand_epsilon(j+1,ih,it+1), ne, ie_sim(j+1,ih,it+1))
                            else
                                if ( rand_upsilon(j,ih,it) < upsilon_d(iz) ) then
                                    in_sim(j+1,ih,it+1) = 3
                                    ie_sim(j+1,ih,it+1) = e
                                else
                                    in_sim(j+1,ih,it+1) = 2
                                    ie_sim(j+1,ih,it+1) = e
                                end if
                            end if
                        else
                            if ( rand_unemp(j,ih,it) <= lambda(iz) ) then
                                in_sim(j+1,ih,it+1) = 1
                                call future_index(cum_pi_e(e,:), rand_epsilon(j+1,ih,it+1), ne, ie_sim(j+1,ih,it+1))
                            else
                                in_sim(j+1,ih,it+1) = 3
                                ie_sim(j+1,ih,it+1) = e
                            end if
                        end if
                    end if
                    
                    !!!! Future Credit State !!!!
                    if ( s == 1 ) then
                        if ( dval <= 0.5_rp ) then
                            is_sim(j+1,ih,it+1) = 1
                        else
                            is_sim(j+1,ih,it+1) = 2
                            ia_sim(j+1,ih,it+1) = nd
                        end if
                    else
                        if ( rand_credit(j,ih,it) <= theta ) then
                            is_sim(j+1,ih,it+1) = 1
                        else
                            is_sim(j+1,ih,it+1) = 2
                        end if
                    end if
                        
                    !!!! Future Discount Factor !!!!
                    ib_sim(j+1,ih,it+1) = b
                        
                    !!!! Household Income !!!!
                    if ( in_sim(j,ih,it) == 1 ) then
                        inc_isim(j,ih,it) = (one-tau)*wage*egrid(e)*gamma_j(j) + y_bar
                    elseif ( in_sim(j,ih,it) == 2 .and. h_isim(j,ih,it) > 0.50_rp ) then
                        if ( upsilon_r(iz)*wage*egrid(e)*gamma_j(j) >= upsilon_bar(iz) ) then
                            inc_isim(j,ih,it) = upsilon_bar(iz) + y_bar
                        else
                            inc_isim(j,ih,it) = upsilon_r(iz)*wage*egrid(e)*gamma_j(j) + y_bar
                        end if
                    else
                        inc_isim(j,ih,it) = y_bar
                    end if
                    
                end if
                    
            end if
                
        end do
    end do
    
end do
    
!!!! Results - Individual Simulation !!!!
it = 0
count_rr = zero
count_quarter = zero; dc_q(:) = zero; dy_q(:) = zero; balance_q(:) = zero; bank_q(:) = zero
dc_mean = zero; count_mean = zero; return_i = zero
            
do t = 980,1188
    it = it + 1
    iz = iz_sim(t)
        
    !!!! Quarterly Statistics after Job Separation !!!!
    if ( t >= 1160 .and. t <= 1170 ) then
        do j = 2,(nj-7)
            do ih = 1,nh
                    
                if ( in_sim(j-1,ih,it-1) == 1 .and. in_sim(j,ih,it) == 2 .and. h_isim(j,ih,it) > 0.5_rp ) then
                    count_quarter = count_quarter + one
                    credit_minus = maxval( (/ -q_isim(j-1,ih,it-1)*g_isim(j-1,ih,it-1), zero /) )
                        
                    return_i = zero
                    do j2 = 1,8
                        dc_q(j2) = dc_q(j2) + ( c_isim(j+j2-2,ih,it+j2-2) - c_isim(j-1,ih,it-1) )/c_isim(j-1,ih,it-1)
                        dy_q(j2) = dy_q(j2) + ( inc_isim(j+j2-2,ih,it+j2-2) - inc_isim(j-1,ih,it-1) )/inc_isim(j-1,ih,it-1)
                        credit_plus = maxval( (/ -q_isim(j+j2-2,ih,it+j2-2)*g_isim(j+j2-2,ih,it+j2-2), zero /) )
                        balance_q(j2) = balance_q(j2) + credit_plus
                            
                        if ( d_isim(j+j2-2,ih,it+j2-2) > 0.5_rp ) then
                            bank_q(j2) = bank_q(j2) + one
                        end if
                            
                        if ( in_sim(j+j2-2,ih,it+j2-2) == 1 .and. j2 > 1 ) then
                            return_i = one
                        end if
                            
                        if ( in_sim(j+j2-2,ih,it+j2-2) > 1 .and. h_isim(j+j2-2,ih,it+j2-2) > 0.5_rp .and. return_i < 0.5_rp ) then
                            count_mean = count_mean + one
                            dc_mean = dc_mean + c_isim(j+j2-2,ih,it+j2-2)/c_isim(j-1,ih,it-1)
                        end if
                    end do
                end if
                    
            end do
        end do
    end if
        
    !!!! Additional Statistics !!!!
    if ( t >= 1160 .and. t <= 1170 ) then
        do j = 2,(nj-3)
            do ih = 1,nh
                    
                if ( in_sim(j-1,ih,it-1) == 1 .and. in_sim(j,ih,it) == 2 .and. h_isim(j,ih,it) > 0.5_rp ) then        
                    open(unit=21, file='Individual_Simulation.txt')
                    write(unit=21, fmt='(f10.6, f10.6)') (c_isim(j,ih,it) - c_isim(j-1,ih,it-1))/c_isim(j-1,ih,it-1), (inc_isim(j,ih,it) - inc_isim(j-1,ih,it-1))/inc_isim(j-1,ih,it-1)
                end if
                        
                !!!! Change in Credit first Quarter of Unemployment !!!!
                if ( in_sim(j-1,ih,it-1) == 1 .and. in_sim(j,ih,it) == 2 .and. h_isim(j,ih,it) > 0.5_rp .and. is_sim(j-1,ih,it-1) == 1 ) then
                    earnings_minus = inc_isim(j-1,ih,it-1)
                    earnings_plus = inc_isim(j,ih,it)
                    debt_minus = maxval( (/ -q_isim(j-1,ih,it-1)*g_isim(j-1,ih,it-1), zero /) )
                    debt_plus = maxval( (/ -q_isim(j,ih,it)*g_isim(j,ih,it), zero /) )
                    if ( earnings_plus < earnings_minus ) then
                        count_rr = count_rr + one
                        if ( debt_plus > debt_minus ) then
                            count_inc = count_inc + one
                        elseif (debt_plus < debt_minus ) then
                            count_dec = count_dec + one
                        end if          
                    end if
                end if
                    
            end do
        end do
    end if
        
end do
    
write(unit=22, fmt=*), ""
write(unit=22, fmt=*), "-------------------------------------------------------------------"
write(unit=22, fmt=*), "   Variables after Job Separation   "
write(unit=22, fmt=*), ""
write(unit=22, fmt='(1x, a6, 5x, a5, 7x, a6, 5x, a7, 5x, a10)'), "Period", "Cons.", "Income", "Balance", "Bankruptcy"
write(unit=22, fmt='(1x, a6, 3x, f8.4, 4x, f8.4, 4x, f8.4, 4x, f8.4)'), "t=-1:", dc_q(1)/count_quarter, dy_q(1)/count_quarter, balance_q(1)/count_quarter, bank_q(1)/count_quarter
write(unit=22, fmt='(1x, a6, 3x, f8.4, 4x, f8.4, 4x, f8.4, 4x, f8.4)'), "t= 0:", dc_q(2)/count_quarter, dy_q(2)/count_quarter, balance_q(2)/count_quarter, bank_q(2)/count_quarter
write(unit=22, fmt='(1x, a6, 3x, f8.4, 4x, f8.4, 4x, f8.4, 4x, f8.4)'), "t= 1:", dc_q(3)/count_quarter, dy_q(3)/count_quarter, balance_q(3)/count_quarter, bank_q(3)/count_quarter
write(unit=22, fmt='(1x, a6, 3x, f8.4, 4x, f8.4, 4x, f8.4, 4x, f8.4)'), "t= 2:", dc_q(4)/count_quarter, dy_q(4)/count_quarter, balance_q(4)/count_quarter, bank_q(4)/count_quarter
write(unit=22, fmt='(1x, a6, 3x, f8.4, 4x, f8.4, 4x, f8.4, 4x, f8.4)'), "t= 3:", dc_q(5)/count_quarter, dy_q(5)/count_quarter, balance_q(5)/count_quarter, bank_q(5)/count_quarter
write(unit=22, fmt='(1x, a6, 3x, f8.4, 4x, f8.4, 4x, f8.4, 4x, f8.4)'), "t= 4:", dc_q(6)/count_quarter, dy_q(6)/count_quarter, balance_q(6)/count_quarter, bank_q(6)/count_quarter
write(unit=22, fmt='(1x, a6, 3x, f8.4, 4x, f8.4, 4x, f8.4, 4x, f8.4)'), "t= 5:", dc_q(7)/count_quarter, dy_q(7)/count_quarter, balance_q(7)/count_quarter, bank_q(7)/count_quarter
write(unit=22, fmt='(1x, a6, 3x, f8.4, 4x, f8.4, 4x, f8.4, 4x, f8.4)'), "t= 6:", dc_q(8)/count_quarter, dy_q(8)/count_quarter, balance_q(8)/count_quarter, bank_q(8)/count_quarter
write(unit=22, fmt=*), "avg cons= ", dc_mean/count_mean
    
write(unit=22, fmt=*), ""
write(unit=22, fmt=*), "-------------------------------------------------------------------"
write(unit=22, fmt=*), "  Credit Change in Quarter of Unemployment  "
write(unit=22, fmt=*), ""
write(unit=22, fmt='(1x, a8, 4x, a8, 3x, a3)'), "Increase", "Decrease", "n_t"
write(unit=22, fmt='(f8.4, 4x, f8.4, 4x, f8.1)'), count_inc/count_rr, count_dec/count_rr, count_rr
    
write(unit=22, fmt=*), ""
write(unit=22, fmt=*), "-------------------------------------------------------------------"
write(unit=22, fmt=*), "  Average MPC  "
write(unit=22, fmt=*), ""
write(unit=22, fmt='(3x, a3)'), "MPC"
write(unit=22, fmt='(1x, f7.4, 4x, f7.4, 4x, f7.4)'), mpc_mean/count_mpc
    
!! Individual Simulation
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!Display Stationary Results - Good Aggregate State
    
!!!! HP Filter Matrix !!!!
HP(:,:) = zero
HP_inv(:,:) = zero
HP(1,1) = one + zeta;           HP(nt-300,nt-300) = one + zeta
HP(1,2) = -2.0_rp*zeta;         HP(nt-300,nt-301) = -2.0_rp*zeta
HP(1,3) = zeta;                 HP(nt-300,nt-302) = zeta
HP(2,1) = -2.0_rp*zeta;         HP(nt-301,nt-300) = -2.0_rp*zeta
HP(2,2) = one + 5.0_rp*zeta;    HP(nt-301,nt-301) = one + 5.0_rp*zeta
HP(2,3) = -4.0_rp*zeta;         HP(nt-301,nt-302) = -4.0_rp*zeta
HP(2,4) = zeta;                 HP(nt-301,nt-303) = zeta;
do t = 3,nt-302
    HP(t,t-2) = zeta
    HP(t,t-1) = -4.0_rp*zeta
    HP(t,t) = one + 6.0_rp*zeta
    HP(t,t+1) = -4.0_rp*zeta
    HP(t,t+2) = zeta
end do
call matrix_inversion(nt-300, HP, HP_inv)
    
!!!! HP Filtered Data !!!!
call matrix_multiplication(nt-300, nt-300, 1, HP_inv, log(lambda_sim(301:nt)), hp_lambda)
call matrix_multiplication(nt-300, nt-300, 1, HP_inv, log(xi_sim(301:nt)), hp_xi)
call matrix_multiplication(nt-300, nt-300, 1, HP_inv, log(Y_sim(301:nt)), hp_Y)
call matrix_multiplication(nt-300, nt-300, 1, HP_inv, log(C_sim(301:nt)), hp_C)
call matrix_multiplication(nt-300, nt-300, 1, HP_inv, log(I_sim(301:nt)), hp_I)
call matrix_multiplication(nt-300, nt-300, 1, HP_inv, log(-D_sim(301:nt)), hp_D)
call matrix_multiplication(nt-300, nt-300, 1, HP_inv, log(B_sim(301:nt)), hp_B)
call matrix_multiplication(nt-300, nt-300, 1, HP_inv, log(M_sim(301:nt)), hp_M)
call matrix_multiplication(nt-300, nt-300, 1, HP_inv, log(U_sim(301:nt)), hp_U)
call matrix_multiplication(nt-300, nt-300, 1, HP_inv, log(P_sim(301:nt)), hp_P)
do t = 1,(nt-300)
    hp_lambda(t) = log(lambda_sim(t+300)) - hp_lambda(t)
    hp_xi(t) = log(xi_sim(t+300)) - hp_xi(t)
    hp_Y(t) = log(Y_sim(t+300)) - hp_Y(t)
    hp_C(t) = log(C_sim(t+300)) - hp_C(t)
    hp_I(t) = log(I_sim(t+300)) - hp_I(t)
    hp_D(t) = log(-D_sim(t+300)) - hp_D(t)
    hp_B(t) = log(B_sim(t+300)) - hp_B(t)
    hp_M(t) = log(M_sim(t+300)) - hp_M(t)
    hp_U(t) = log(U_sim(t+300)) - hp_U(t)
    hp_P(t) = log(P_sim(t+300)) - hp_P(t)
end do
    
!!!! First Moments !!!!
mean_Z = sum(Z_sim(301:nt))/dble(nt-300)
mean_lambda = sum(lambda_sim(301:nt))/dble(nt-300)
mean_xi = sum(xi_sim(301:nt))/dble(nt-300)
mean_K = sum(K_sim(301:nt))/dble(nt-300)
mean_L = sum(L_sim(301:nt))/dble(nt-300)
mean_Y = sum(Y_sim(301:nt))/dble(nt-300)
mean_C = sum(C_sim(301:nt))/dble(nt-300)
mean_I = sum(I_sim(301:nt))/dble(nt-300)
mean_G = sum(G_sim(301:nt))/dble(nt-300)
mean_D = sum(D_sim(301:nt))/dble(nt-300)
mean_B = sum(B_sim(301:nt))/dble(nt-300)
mean_U = sum(U_sim(301:nt))/dble(nt-300)
mean_P = sum(P_sim(301:nt))/dble(nt-300)
mean_r = sum(r_sim(301:nt))/dble(nt-300)
mean_w = sum(w_sim(301:nt))/dble(nt-300)
mean_M = sum(M_sim(301:nt))/dble(nt-300)
mean_abar = sum(abar_sim(301:nt))/dble(nt-300)
mean_ebar = sum(ebar_sim(301:nt))/dble(nt-300)
    
!!!! Second Moments !!!!
call standard_deviation(nt-300, sum(hp_Z(:))/dble(nt-300), hp_Z(:), sd_Z)
call standard_deviation(nt-300, sum(hp_lambda(:))/dble(nt-300), hp_lambda(:), sd_lambda)
call standard_deviation(nt-300, sum(hp_xi(:))/dble(nt-300), hp_xi(:), sd_xi)
call standard_deviation(nt-300, sum(hp_Y(:))/dble(nt-300), hp_Y(:), sd_Y)
call standard_deviation(nt-300, sum(hp_C(:))/dble(nt-300), hp_C(:), sd_C)
call standard_deviation(nt-300, sum(hp_I(:))/dble(nt-300), hp_I(:), sd_I)
call standard_deviation(nt-300, sum(hp_D(:))/dble(nt-300), hp_D(:), sd_D)
call standard_deviation(nt-300, sum(hp_B(:))/dble(nt-300), hp_B(:), sd_B)
call standard_deviation(nt-300, sum(hp_U(:))/dble(nt-300), hp_U(:), sd_U)
call standard_deviation(nt-300, sum(hp_P(:))/dble(nt-300), hp_P(:), sd_P)
call standard_deviation(nt-300, sum(hp_M(:))/dble(nt-300), hp_M(:), sd_M)

!!!! Cross-Correlation !!!!
call correlation_coefficient(nt-300, sum(hp_Y(:))/dble(nt-300), sum(hp_Y(:))/dble(nt-300), HP_Y(:), HP_Y(:), corr_Y)
call correlation_coefficient(nt-300, sum(hp_C(:))/dble(nt-300), sum(hp_Y(:))/dble(nt-300), HP_C(:), HP_Y(:), corr_C)
call correlation_coefficient(nt-300, sum(hp_I(:))/dble(nt-300), sum(hp_Y(:))/dble(nt-300), HP_I(:), HP_Y(:), corr_I)
call correlation_coefficient(nt-300, sum(hp_D(:))/dble(nt-300), sum(hp_Y(:))/dble(nt-300), HP_D(:), HP_Y(:), corr_D)
call correlation_coefficient(nt-300, sum(hp_B(:))/dble(nt-300), sum(hp_Y(:))/dble(nt-300), HP_B(:), HP_Y(:), corr_B)
call correlation_coefficient(nt-300, sum(hp_U(:))/dble(nt-300), sum(hp_Y(:))/dble(nt-300), HP_U(:), HP_Y(:), corr_U)
call correlation_coefficient(nt-300, sum(hp_P(:))/dble(nt-300), sum(hp_Y(:))/dble(nt-300), HP_P(:), HP_Y(:), corr_P)
call correlation_coefficient(nt-300, sum(hp_M(:))/dble(nt-300), sum(hp_Y(:))/dble(nt-300), HP_M(:), HP_Y(:), corr_M)

!!!! Aggregate Variables !!!!
if ( switch_results == 1 ) then
    write(unit=22, fmt=*), ""
    write(unit=22, fmt=*), ""
    write(unit=22, fmt=*), "----------------------------------------------------------------------"
    write(unit=22, fmt=*), "----------------- Exogenous Stochastic Components --------------------"
    write(unit=22, fmt=*), ""
    write(unit=22, fmt='(24x, a4, 8x, a7)'), "Mean", "Std Dev"
    write(unit=22, fmt='(1x, a5, 15x, f8.4, 5x, a4)'), "TFP: ", mean_Z, "  NA"
    write(unit=22, fmt='(1x, a13, 7x, f8.4, 5x, f8.4)'), "Job Finding: ", mean_lambda, sd_lambda
    write(unit=22, fmt='(1x, a17, 3x, f8.4, 5x, f8.4)'), "Job Separations: ", mean_xi, sd_xi
end if

!!!! Cyclical Moments !!!!
if ( switch_results == 1 ) then
    write(unit=22, fmt=*), ""
    write(unit=22, fmt=*), ""
    write(unit=22, fmt=*), "----------------------------------------------------------------------"
    write(unit=22, fmt=*), "------------------------- Cyclical Moments ---------------------------"
    write(unit=22, fmt=*), ""
    write(unit=22, fmt='(24x, a4, 8x, a7, 4x, a10, 4x, a9)'), "Mean", "Std Dev", "Cross-Corr"
    write(unit=22, fmt='(1x, a8, 12x, f8.4, 5x, f8.4, 5x, f8.4)'), "Output: ", mean_Y, sd_Y, corr_Y
    write(unit=22, fmt='(1x, a13, 7x, f8.4, 5x, f8.4, 5x, f8.4)'), "Consumption: ", mean_C, sd_C, corr_C
    write(unit=22, fmt='(1x, a12, 8x, f8.4, 5x, f8.4, 5x, f8.4)'), "Investment: ", mean_I, sd_I, corr_I
    write(unit=22, fmt='(1x, a12, 8x, f8.4, 5x, f8.4, 5x, f8.4)'), "Uns Credit: ", mean_D, sd_D, corr_D
    write(unit=22, fmt='(1x, a14, 6x, f8.4, 5x, f8.4, 5x, f8.4)'), "Bankruptcies: ", mean_B, sd_B, corr_B
    write(unit=22, fmt='(1x, a16, 4x, f8.4, 5x, f8.4, 5x, f8.4)'), "Credit Premium: ", mean_M, sd_M, corr_M
    write(unit=22, fmt='(1x, a14, 6x, f8.4, 5x, f8.4, 5x, f8.4)'), "Unemployment: ", mean_U, sd_U, corr_U
    write(unit=22, fmt='(1x, a15, 5x, f8.4, 5x, f8.4, 5x, f8.4)'), "Participation: ", mean_P, sd_P, corr_P
end if

!!!! Calibration Targets !!!!
if ( switch_results == 1 ) then
    write(unit=22, fmt=*), ""
    write(unit=22, fmt=*), ""
    write(unit=22, fmt=*), "----------------------------------------------------------------------"
	write(unit=22, fmt=*), "------------- Calibration Targets (Annualized 1980-2019) -------------"
	write(unit=22, fmt=*), ""
	write(unit=22, fmt='(a9, 3x, a5, 6x, a6, 19x, a4, 7x, a5 )'), "Parameter", "Value", "Target", "Data", "Model"
	write(unit=22, fmt='(1x, a6, 3x, f7.3, 4x, a17, 8x, f7.3, 5x, f7.3 )'), "beta_h", beta(1), "Capital-GDP Ratio", 224.07, (mean_K/(4.0_rp*mean_Y))*100.0_rp
    write(unit=22, fmt='(1x, a5, 4x, f7.3, 4x, a17, 8x, f7.3, 5x, f7.3 )'), "delta", delta, "Inv-Capital Ratio", 7.89, ((4.0_rp*mean_I)/mean_K)*100.0_rp
	write(unit=22, fmt='(1x, a6, 3x, f7.3, 4x, a14, 11x, f7.3, 5x, f7.3 )'), "beta_l", beta(2), "Debt-GDP Ratio", 4.83, -(mean_D/(4.0_rp*mean_Y))*100.0_rp
    write(unit=22, fmt='(1x, a7, 2x, f7.3, 4x, a17, 8x, f7.3, 5x, f7.3 )'), "pi_beta", pi_beta, "Share of LF Borrowing", 41.37, 100.0_rp*Dpop_P_mean/dble(nt-300)
	write(unit=22, fmt='(1x, a5, 4x, f7.3, 4x, a22, 3x, f7.3, 5x, f7.3 )'), "chi_b", chi_b, "Chapter 7 Bankruptcies", 0.88, 4.0_rp*mean_B*100.0_rp
    write(unit=22, fmt='(1x, a4, 5x, f7.3, 4x, a14, 11x, f7.3, 5x, f7.3 )'), "iota", iota_g, "Credit Premium", 11.73, mean_M
    write(unit=22, fmt='(1x, a5, 4x, f7.3, 4x, a13, 12x, f7.3, 5x, f7.3 )'), "y_bar", y_bar, "Gov-GDP Ratio", 5.38, (mean_G/mean_Y)*100.0_rp
    write(unit=22, fmt='(1x, a5, 4x, f7.3, 4x, a18, 7x, f7.3, 5x, f7.3 )'), "chi_s", chi_s, "Participation Rate", 73.73, mean_P
    write(unit=22, fmt='(1x, a8, 1x, f7.3, 4x, a19, 6x, f7.3, 5x, f7.3 )'), "lambda_g", lambda_g, "Mean Unemp Duration", 20.93, 13.0_rp/mean_lambda
	write(unit=22, fmt='(1x, a4, 5x, f7.3, 4x, a17, 8x, f7.3, 5x, f7.3 )'), "xi_g", xi_g, "Unemployment Rate", 6.32, mean_U
    write(unit=22, fmt='(1x, a4, 5x, f7.3, 4x, a13, 12x, f7.3, 5x, f7.3 )'), "xi_b", xi_b, "SD Unemp Rate", 11.12, 100.0_rp*sd_U
    write(unit=22, fmt='(1x, a8, 1x, f7.3, 4x, a15, 10x, f7.3, 5x, f7.3 )'), "lambda_b", lambda_b, "SD_lambda/SD_xi", 1.62, sd_lambda/sd_xi
    write(unit=22, fmt='(1x, a3, 6x, f7.3, 4x, a6, 19x, f7.3, 5x, f7.3 )'), "z_b", z_b, "SD GDP", 1.23, 100.0_rp*sd_Y
    write(unit=22, fmt='(1x, a8, 1x, f7.3, 4x, a16, 9x, f7.3, 5x, f7.3)'), "avg_earn", e_bar, "Average Earnings", e_bar, mean_ebar
end if


end program BCI_Main